code refactoring.

checked if the theoretical fundamental diagram is correct
(yes it is!)

asep.py

 import numpy as np
 import matplotlib.pyplot as plt
 from mpl_toolkits.axes_grid1 import make_axes_locatable
-import time, random 
-import logging, types, argparse
+import time
+import logging
+import argparse
 
 logfile = 'log.dat'
 logging.basicConfig(filename=logfile, level=logging.INFO, format='%(asctime)s - %(levelname)s - %(message)s')
-def getParserArgs():
+
+
+def get_parser_args():
     parser = argparse.ArgumentParser(description='ASEP - TASEP')
     parser.add_argument("-n", "--np", type=int, default=10, help='number of agents (default 10)')
     parser.add_argument("-N", "--nr", type=int, default=1, help='number of runs (default 1)')
+    parser.add_argument("-m", "--ms", type=int, default=100, help='max simulation steps (default 100)')
     parser.add_argument("-p", "--plotP", action='store_const', const=1, default=0, help='plot Pedestrians')
     parser.add_argument("-r", "--shuffle", action='store_const', const=1, default=0, help='random shuffle')
     parser.add_argument("-v", "--reverse", action='store_const', const=1, default=0, help='reverse sequential update')
-    parser.add_argument("-l", "--log" , type=argparse.FileType('w'), default='log.dat', help="log file (default log.dat)")
+    parser.add_argument("-l", "--log", type=argparse.FileType('w'), default='log.dat',
+                        help="log file (default log.dat)")
     args = parser.parse_args()
     return args
 
 
-def init_peds(N, n_cells):
+def init_cells(num_peds, num_cells):
     """
     distribute N pedestrians in box 
     """
-    if N>n_cells: 
-        N = n_cells
-    PEDS = np.ones(N, int) # pedestrians    
-    EMPTY_CELLS_in_BOX = np.zeros( n_cells - N, int) # the rest of cells in the box
-    PEDS = np.hstack((PEDS,  EMPTY_CELLS_in_BOX)) # put 0s and 1s together
-    np.random.shuffle(PEDS) # shuffle them
-    return PEDS
+    if num_peds > num_cells:
+        num_peds = num_cells
+        
+    cells = np.ones(num_peds, int)  # pedestrians
+    zero_cells = np.zeros(num_cells - num_peds, int)  # the rest of cells in the box
+    cells = np.hstack((cells, zero_cells))  # put 0s and 1s together
+    np.random.shuffle(cells)  # shuffle them
+    return cells
 
 
-def plot_peds(peds, walls, i):
+def plot_cells(state_cells, walls_inf, i):
     """
-    plot the cells. we need to make 'bad' walls to better visualize the cells
+    plot the actual state of the cells. we need to make 'bad' walls to better visualize the cells
     """
-    walls = walls * np.Inf
-    P = np.vstack((walls,peds))
-    P = np.vstack((P,walls))
+    walls_inf = walls_inf * np.Inf
+    tmp_cells = np.vstack((walls_inf, state_cells))
+    tmp_cells = np.vstack((tmp_cells, walls_inf))
     fig = plt.figure()
     ax = fig.add_subplot(111)
     ax.cla()
     cmap = plt.get_cmap("gray")
     cmap.set_bad(color='k', alpha=0.8)
-    im = ax.imshow(P, cmap=cmap, vmin=0, vmax=1,interpolation = 'nearest')
+    im = ax.imshow(tmp_cells, cmap=cmap, vmin=0, vmax=1, interpolation='nearest')
     divider = make_axes_locatable(ax)
     cax = divider.append_axes("right", size="1%", pad=0.1)
     plt.colorbar(im, cax=cax, ticks=[0, 1])
     ax.set_axis_off()
-    N = sum(peds)
-    S = "t: %3.3d | n: %d\n"%(i,N)
-    plt.title("%20s"%S,rotation=0,fontsize=10,verticalalignment='bottom')
+    N = sum(state_cells)
+    text = "t: %3.3d | n: %d\n" % (i, N)
+    plt.title("%20s" % text, rotation=0, fontsize=10, verticalalignment='bottom')
     plt.savefig("figs/peds%.5d.png" % i, dpi=100, facecolor='lightgray')
 
-def print_logs(N_pedestrians, width,t, dt, nruns, Dt, vel, d):
+
+def print_logs(num_pedestrians,  system_width, simulation_steps, evac_time, total_runtime, nruns, vel, d):
     """
-    print some infos to the screen
+    print some information to the screen
+    :rtype : none
     """
-    print ("Simulation space (%.2f x 1) m^2"%(width))
-    print ("Mean Evacuation time: %.2f s, runs: %d"%((t*dt)/nruns, nruns))
-    print ("Total Run time: %.2f s"%(Dt))
-    print ("Factor: %.2f s"%( dt*t/Dt ) )
-    print("--------------------------")
-    print ("N %d   mean_velocity  %.2f [m/s]   density  %.2f [1/m]"%(N_pedestrians, vel, d))
-    print("--------------------------")
-    
-def asep_parallel(cells):
-    neighbors = np.roll(cells,-1)
-    dim = len(cells)
-    nmove = 0
+    print('\n')
+    print ('Simulation space (%.2f x 1) m^2' % system_width)
+    print ('Mean Evacuation time: %.2f s, runs: %d' % ((evac_time * dt) / nruns, nruns))
+    print ('max simulation steps %d' % simulation_steps)
+    print ('Total Run time: %.2f s' % total_runtime)
+    print ('Factor: %.2f s' % (dt * evac_time / total_runtime))
+    print('--------------------------')
+    print ('N %d   mean_velocity  %.2f [m/s]   density  %.2f [1/m]' % (num_pedestrians, vel, d))
+    print('--------------------------')
+
+
+def asep_parallel(actual_cells):
+    neighbors = np.roll(actual_cells, -1)
+    assert isinstance(actual_cells, np.ndarray)
+    dim = len(actual_cells)
+    num_move = 0
     tmp_cells = np.zeros(dim)
-    for i,j in  enumerate(cells):
+    for i, j in enumerate(actual_cells):
         if j and not neighbors[i]:
             tmp_cells[i], tmp_cells[(i + 1) % dim] = 0, 1
-            nmove += 1
+            num_move += 1
         elif j:
             tmp_cells[i] = 1
-    return tmp_cells, nmove
+    return tmp_cells, num_move
+
 
 if __name__ == "__main__":
-    args = getParserArgs() # get arguments
+    args = get_parser_args()  # get arguments
     # init parameters
     N_pedestrians = args.np
     shuffle = args.shuffle
     reverse = args.reverse
     drawP = args.plotP
     #######################################################
-    MAX_STEPS = 20 # simulation time
-    nruns = args.nr  # repeate simulations, for TASEP 
-    steps = range(MAX_STEPS)
+    max_steps = args.ms  # simulation time
+    num_runs = args.nr  # repeat simulations, for TASEP
+    steps = range(max_steps)
     cellSize = 0.4  # m
-    vmax= 1.2 # m/s
-    dt = cellSize/vmax  # time step
-    width = 10 # m
+    max_velocity = 1.2  # m/s
+    dt = cellSize / max_velocity  # time step
+    width = 40  # m
     n_cells = int(width / cellSize)  # number of cells
     if N_pedestrians >= n_cells:
         N_pedestrians = n_cells - 1
-        print("warning: maximum of %d cells are allowed"%N_pedestrians)
+        print('warning: maximum of %d cells are allowed' % N_pedestrians)
     else:
-        print("info: n_pedestrians=%d (max_pedestrians=%d)"%(N_pedestrians,n_cells))
+        print('info: n_pedestrians=%d (max_pedestrians=%d)' % (N_pedestrians, n_cells))
     #######################################################
     t1 = time.time()
-    tsim = 0
+    simulation_time = 0
     density = float(N_pedestrians) / width
     walls = np.ones(n_cells)
     velocities = []  # velocities over all runs
-    for n in range(nruns):
-        peds = init_peds(N_pedestrians, n_cells)
+    for n in range(num_runs):
+        cells = init_cells(N_pedestrians, n_cells)
         velocity = 0
-        for t in steps: # simulation loop
+        for step in steps:  # simulation loop
             if drawP:
-                plot_peds(peds, walls, t)
+                plot_cells(cells, walls, step)
 
-            peds, nmove = asep_parallel(peds)
-            v = nmove*vmax/float(N_pedestrians)
+            cells, num_moves = asep_parallel(cells)
+            v = num_moves * max_velocity / float(N_pedestrians)
             velocity += v
 
-        velocity /= MAX_STEPS
+        velocity /= max_steps
         velocities.append(velocity)
-        tsim += t
-        print ("\trun: %3d ----  vel: %.2f |  den: %.2f"%(n, velocity,  density))
+        simulation_time += max_steps
+        print ('\t run: %3d ----  vel: %.2f |  den: %.2f' % (n, velocity, density))
     t2 = time.time()
     mean_velocity = np.mean(velocities)
-    print_logs(N_pedestrians, width, tsim, dt, nruns, t2-t1, mean_velocity, density)
-    
-    
+    print_logs(N_pedestrians, width, max_steps, simulation_time, t2-t1, num_runs, mean_velocity, density)

make_fd.py

@@ -3,29 +3,31 @@ import subprocess
 import matplotlib.pyplot as plt
 
 # ----------------------------------------
-nruns = 30
-maxpeds = 30
-npeds = range(1, maxpeds)
-stdout = open("stdout.txt","w")
+num_runs = 30
+max_pedestrians = 120
+sim_steps = 1000
+pedestrians = range(1, max_pedestrians)
+filename = open("stdout.txt", "w")
 # ----------------------------------------
-for n in npeds:
-    print("run asep.py with -n %3.3d -N %3.3d"%(n,nruns))
-    subprocess.call(["python", "asep.py", "-n" "%d"%n, "-N", "%d"%nruns], stdout=stdout)
+for n in pedestrians:
+    print("run asep.py with -n %3.3d -N %3.3d -m % 3.4d" % (n, num_runs, sim_steps))
+    subprocess.call(["python", "asep.py", "-n" "%d" % n, "-N", "%d" % num_runs, "-m", "%d" % sim_steps],
+                    stdout=filename)
 # ----------------------------------------
-stdout.close()
+filename.close()
 
 velocities = []
 densities = []
 # the line should be something like this
 # N 1   mean_velocity  1.20 [m/s]   density  0.10 [1/m]
-stdout = open("stdout.txt","r")
-for line in stdout:
+filename = open("stdout.txt", "r")
+for line in filename:
     if line.startswith("N"):
         line = line.split()
         velocities.append(float(line[3]))
         densities.append(float(line[6]))
 
-stdout.close()
+filename.close()
 # -------- plot FD ----------
 # rho vs v
 fig = plt.figure()
@@ -44,6 +46,6 @@ plt.ylabel(r"$J\, [s^{-1}]$", size=20)
 fig.tight_layout()
 print("\n")
 for end in ["pdf", "png", "eps"]:
-    figname = "asep_fd.%s"%end
-    print("result in %s"%figname)
-    plt.savefig(figname)
+    figure_name = "asep_fd.%s" % end
+    print("result written in %s" % figure_name)
+    plt.savefig(figure_name)
\ No newline at end of file