Skip to main content
`Ctrl`+`K`

# IPF Performance

The use of iterative proportional fitting (IPF) is quite common on processes involving doubly-constraining matrices, such as synthetic gravity models and fractional split models (aggregate destination-choice models).

As this is a commonly used algorithm, we have implemented it in Cython, where we can take full advantage of multi-core CPUs. We have also implemented the ability of using both 32-bit and 64-bit floating-point seed matrices, which has direct impact on cache use and consequently computational performance.

In this section, we compare the runtime of AequilibraE’s current implementation of IPF, with a general IPF algorithm written in pure Python, available here.

The figure below compares AequilibraE’s IPF runtime with one core with the benchmark Python code. From the figure below, we can notice that the runtimes were practically the same for the instances with 1,000 zones or less. As the number of zones increases, AequilibraE demonstrated to be slightly faster than the benchmark python code, while applying IPF to a 32-bit NumPy array (np.float32) was significantly faster.

It’s worth mentioning that the user can set up a threshold for AequilibraE’s IPF function, as well as use more than one core to speed up the fitting process.

AequilibraE’s IPF runtime

As IPF is an embarrassingly-parallel workload, it is more relevant to look at the performance of the AequilibraE implementations, starting by comparing the implementation performance for inputs in 32 vs 64 bits using 32 threads.

AequilibraE’s IPF runtime 32 vs 64 bits

The difference is staggering, with the 32-bit implementation being twice as fast as the 64-bit one for large matrices. It is also worth noting that differences in results between the outputs between these two versions are incredibly small (RMSE < 1.1e-10), and therefore unlikely to be relevant in most applications.

We can also look at performance gain across matrix sizes and number of cores, and it becomes clear that the 32-bit version scales significantly better than its 64-bit counterpart, showing significant performance gains up to 16 threads, while the latter stops showing much improvement beyond 8 threads, likely due to limitations on cache size.

number of cores used in IPF for 64 bit matrices

number of cores used in IPF for 32 bit matrices

In conclusion, AequilibraE’s IPF implementation is over 11 times faster than its pure Python counterpart for large matrices on a workstation, largely due to the use of Cython and multi-threading, but also due to the use of a 32-bit version of the algorithm.

These tests were run on a Threadripper 3970x (released in 2019) workstation with 32 cores (64 threads) @ 3.7 GHz and 256 Gb of RAM. The code is provided below for reference.

## Reference code

    from copy import deepcopy
    from time import perf_counter
    import numpy as np
    import pandas as pd
    from aequilibrae.distribution.ipf_core import ipf_core
    from tqdm import tqdm

    # From:
    # https://github.com/joshchea/python-tdm/blob/master/scripts/CalcDistribution.py

    def CalcFratar(ProdA, AttrA, Trips1, maxIter=10):
        '''Calculates fratar trip distribution
           ProdA = Production target as array
           AttrA = Attraction target as array
           Trips1 = Seed trip table for fratar
           maxIter (optional) = maximum iterations, default is 10
           Returns fratared trip table
        '''
        # print('Checking production, attraction balancing:')
        sumP = ProdA.sum()
        sumA = AttrA.sum()
        # print('Production: ', sumP)
        # print('Attraction: ', sumA)
        if sumP != sumA:
            # print('Productions and attractions do not balance, attractions will be scaled to productions!')
            AttrA = AttrA*(sumP/sumA)
        else:
            pass
            # print('Production, attraction balancing OK.')
        # Run 2D balancing --->
        for balIter in range(0, maxIter):
            ComputedProductions = Trips1.sum(1)
            ComputedProductions[ComputedProductions == 0] = 1
            OrigFac = (ProdA/ComputedProductions)
            Trips1 = Trips1*OrigFac[:, np.newaxis]

            ComputedAttractions = Trips1.sum(0)
            ComputedAttractions[ComputedAttractions == 0] = 1
            DestFac = (AttrA/ComputedAttractions)
            Trips1 = Trips1*DestFac
        return Trips1

    mat_sizes = [500, 750, 1000, 1500, 2500, 5000, 7500, 10000, 15000]

    #Benchmarking
    bench_data = []
    cores = 1
    repetitions = 5
    iterations = 100
    for zones in mat_sizes:
        for repeat in tqdm(range(repetitions), f"Repetitions for zone size {zones}"):
            mat1 = np.random.rand(zones, zones)
            target_prod = np.random.rand(zones)
            target_atra = np.random.rand(zones)
            target_atra *= target_prod.sum()/target_atra.sum()

            aeq_mat = deepcopy(mat1)
            # We use a nonsensical negative tolerance to force it to run all iterations
            # and set warning for non-convergence to false, as we know it won't converge
            t = perf_counter()
            ipf_core(aeq_mat, target_prod, target_atra, max_iterations=iterations, tolerance=-5, cores=cores, warn=False)
            aeqt = perf_counter() - t

            aeq_mat32 = np.array(mat1, np.float32)
            # We now run the same thing with a seed matrix in single-precision (float 32 bits) instead of double as above (64 bits)
            t = perf_counter()
            ipf_core(aeq_mat32, target_prod, target_atra, max_iterations=iterations, tolerance=-5, cores=cores, warn=False)
            aeqt2 = perf_counter() - t

            bc_mat = deepcopy(mat1)
            t = perf_counter()
            x = CalcFratar(target_prod, target_atra, bc_mat, maxIter=iterations)

            bench_data.append([zones, perf_counter() - t, aeqt, aeqt2])

    Repetitions for zone size 500: 100%|██████████| 5/5 [00:01<00:00,  2.60it/s]
    Repetitions for zone size 750: 100%|██████████| 5/5 [00:04<00:00,  1.18it/s]
    Repetitions for zone size 1000: 100%|██████████| 5/5 [00:07<00:00,  1.57s/it]
    Repetitions for zone size 1500: 100%|██████████| 5/5 [00:19<00:00,  3.88s/it]
    Repetitions for zone size 2500: 100%|██████████| 5/5 [00:56<00:00, 11.24s/it]
    Repetitions for zone size 5000: 100%|██████████| 5/5 [03:44<00:00, 44.89s/it]
    Repetitions for zone size 7500: 100%|██████████| 5/5 [08:09<00:00, 97.89s/it]
    Repetitions for zone size 10000: 100%|██████████| 5/5 [14:11<00:00, 170.34s/it]
    Repetitions for zone size 15000: 100%|██████████| 5/5 [32:23<00:00, 388.70s/it]

    bench_df = pd.DataFrame(bench_data, columns=["Zones in the model", "PythonCode", "AequilibraE", "AequilibraE-32bits"])
    bench_df.groupby(["Zones in the model"]).mean().plot.bar()

    <Axes: xlabel='Zones in the model'>

    bench_df.groupby(["Zones in the model"]).mean()

| PythonCode | AequilibraE | AequilibraE-32bits
--- | --- | --- | ---
Zones in the model |  |  |
500 | 0.103166 | 0.140051 | 0.137266
750 | 0.241757 | 0.302931 | 0.297389
1000 | 0.513984 | 0.529727 | 0.516339
1500 | 1.416848 | 1.276844 | 1.156637
2500 | 4.253563 | 3.581645 | 3.333664
5000 | 17.059166 | 14.288760 | 13.257936
7500 | 35.568686 | 32.021563 | 29.666000
10000 | 59.030967 | 57.028204 | 53.155270
15000 | 134.895915 | 130.652199 | 120.535656

    #Benchmarking 32 threads
    bench_data_parallel = []
    cores = 32
    repetitions = 5
    iterations = 100
    for zones in mat_sizes:
        for repeat in tqdm(range(repetitions), f"Repetitions for zone size {zones}"):
            mat1 = np.random.rand(zones, zones)
            target_prod = np.random.rand(zones)
            target_atra = np.random.rand(zones)
            target_atra *= target_prod.sum()/target_atra.sum()

            aeq_mat = deepcopy(mat1)
            # We use a nonsensical negative tolerance to force it to run all iterations
            # and set warning for non-convergence to false, as we know it won't converge
            t = perf_counter()
            ipf_core(aeq_mat, target_prod, target_atra, max_iterations=iterations, tolerance=-5, cores=cores, warn=False)
            aeqt = perf_counter() - t

            aeq_mat32 = np.array(mat1, np.float32)
            # We now run the same thing with a seed matrix in single-precision (float 32 bits) instead of double as above (64 bits)
            t = perf_counter()
            ipf_core(aeq_mat32, target_prod, target_atra, max_iterations=iterations, tolerance=-5, cores=cores, warn=False)
            aeqt2 = perf_counter() - t

            rmse = np.sqrt(np.mean((aeq_mat-aeq_mat32)**2))

            bench_data_parallel.append([zones, aeqt, aeqt2, rmse])

    Repetitions for zone size 500: 100%|██████████| 5/5 [00:01<00:00,  2.70it/s]
    Repetitions for zone size 750: 100%|██████████| 5/5 [00:01<00:00,  2.64it/s]
    Repetitions for zone size 1000: 100%|██████████| 5/5 [00:02<00:00,  2.37it/s]
    Repetitions for zone size 1500: 100%|██████████| 5/5 [00:03<00:00,  1.61it/s]
    Repetitions for zone size 2500: 100%|██████████| 5/5 [00:07<00:00,  1.41s/it]
    Repetitions for zone size 5000: 100%|██████████| 5/5 [00:24<00:00,  4.91s/it]
    Repetitions for zone size 7500: 100%|██████████| 5/5 [00:49<00:00,  9.96s/it]
    Repetitions for zone size 10000: 100%|██████████| 5/5 [01:26<00:00, 17.29s/it]
    Repetitions for zone size 15000: 100%|██████████| 5/5 [03:10<00:00, 38.02s/it]

    bench_df_parallel = pd.DataFrame(bench_data_parallel, columns=["Zones in the model", "AequilibraE", "AequilibraE-32bits", "rmse"])
    bench_df_parallel.groupby(["Zones in the model"]).mean()[[ "AequilibraE", "AequilibraE-32bits"]].plot.bar()

    <Axes: xlabel='Zones in the model'>

    bench_df_parallel.groupby(["Zones in the model"]).mean()

| AequilibraE | AequilibraE-32bits | rmse
--- | --- | --- | ---
Zones in the model |  |  |
500 | 0.170232 | 0.172460 | 1.080659e-10
750 | 0.183474 | 0.183857 | 7.100473e-11
1000 | 0.211237 | 0.191961 | 5.137420e-11
1500 | 0.322709 | 0.252115 | 3.358886e-11
2500 | 0.779715 | 0.522037 | 2.036197e-11
5000 | 2.745004 | 1.760100 | 9.588117e-12
7500 | 5.900431 | 3.170006 | 6.258958e-12
10000 | 10.255305 | 5.456631 | 4.733115e-12
15000 | 22.618665 | 11.860497 | 3.102268e-12

    cores_to_use = [1, 2, 4, 8, 16, 32]

    aeq_data = []
    repetitions = 1
    iterations = 50
    for zones in mat_sizes:
        for cores in tqdm(cores_to_use,f"Zone size: {zones}"):
            for repeat in range(repetitions):
                mat1 = np.random.rand(zones, zones)
                target_prod = np.random.rand(zones)
                target_atra = np.random.rand(zones)
                target_atra *= target_prod.sum()/target_atra.sum()

                aeq_mat = np.array(deepcopy(mat1), np.float32)
                t = perf_counter()
                ipf_core(aeq_mat, target_prod, target_atra, max_iterations=iterations, tolerance=-5, cores=cores, warn=False)
                aeqt = perf_counter() - t

                aeq_data.append([zones, cores, aeqt])

    Zone size: 500: 100%|██████████| 6/6 [00:00<00:00, 12.14it/s]
    Zone size: 750: 100%|██████████| 6/6 [00:00<00:00, 10.20it/s]
    Zone size: 1000: 100%|██████████| 6/6 [00:00<00:00,  6.87it/s]
    Zone size: 1500: 100%|██████████| 6/6 [00:01<00:00,  3.42it/s]
    Zone size: 2500: 100%|██████████| 6/6 [00:04<00:00,  1.32it/s]
    Zone size: 5000: 100%|██████████| 6/6 [00:16<00:00,  2.73s/it]
    Zone size: 7500: 100%|██████████| 6/6 [00:35<00:00,  5.93s/it]
    Zone size: 10000: 100%|██████████| 6/6 [01:02<00:00, 10.46s/it]
    Zone size: 15000: 100%|██████████| 6/6 [02:21<00:00, 23.62s/it]

    aeq_df = pd.DataFrame(aeq_data, columns=["zones", "cores", "time"])
    aeq_df = aeq_df[aeq_df.zones>1000]
    aeq_df = aeq_df.groupby(["zones", "cores"]).mean().reset_index()
    aeq_df = aeq_df.pivot_table(index="zones", columns="cores", values="time")
    for cores in cores_to_use[::-1]:
        aeq_df.loc[:, cores] /= aeq_df[1]
    aeq_df.transpose().plot()
    aeq_df

cores | 1 | 2 | 4 | 8 | 16 | 32
--- | --- | --- | --- | --- | --- | ---
zones |  |  |  |  |  |
1500 | 1.0 | 0.563866 | 0.426711 | 0.270452 | 0.191147 | 0.213856
2500 | 1.0 | 0.527781 | 0.375033 | 0.249776 | 0.156473 | 0.157657
5000 | 1.0 | 0.510183 | 0.278150 | 0.174469 | 0.143175 | 0.134064
7500 | 1.0 | 0.508706 | 0.265835 | 0.157217 | 0.125007 | 0.107014
10000 | 1.0 | 0.507172 | 0.266192 | 0.154844 | 0.120771 | 0.103968
15000 | 1.0 | 0.504696 | 0.259933 | 0.152394 | 0.115052 | 0.101344

    aeq_data = []
    repetitions = 1
    iterations = 50
    for zones in mat_sizes:
        for cores in tqdm(cores_to_use,f"Zone size: {zones}"):
            for repeat in range(repetitions):
                mat1 = np.random.rand(zones, zones)
                target_prod = np.random.rand(zones)
                target_atra = np.random.rand(zones)
                target_atra *= target_prod.sum()/target_atra.sum()

                aeq_mat = np.array(deepcopy(mat1), np.float64)
                t = perf_counter()
                ipf_core(aeq_mat, target_prod, target_atra, max_iterations=iterations, tolerance=-5, cores=cores, warn=False)
                aeqt = perf_counter() - t

                aeq_data.append([zones, cores, aeqt])

    Zone size: 500: 100%|██████████| 6/6 [00:00<00:00, 12.51it/s]
    Zone size: 750: 100%|██████████| 6/6 [00:00<00:00,  9.19it/s]
    Zone size: 1000: 100%|██████████| 6/6 [00:00<00:00,  6.50it/s]
    Zone size: 1500: 100%|██████████| 6/6 [00:01<00:00,  3.07it/s]
    Zone size: 2500: 100%|██████████| 6/6 [00:05<00:00,  1.17it/s]
    Zone size: 5000: 100%|██████████| 6/6 [00:18<00:00,  3.14s/it]
    Zone size: 7500: 100%|██████████| 6/6 [00:42<00:00,  7.10s/it]
    Zone size: 10000: 100%|██████████| 6/6 [01:15<00:00, 12.51s/it]
    Zone size: 15000: 100%|██████████| 6/6 [02:47<00:00, 27.93s/it]

    aeq_df = pd.DataFrame(aeq_data, columns=["zones", "cores", "time"])
    aeq_df = aeq_df[aeq_df.zones>1000]
    aeq_df = aeq_df.groupby(["zones", "cores"]).mean().reset_index()
    aeq_df = aeq_df.pivot_table(index="zones", columns="cores", values="time")
    for cores in cores_to_use[::-1]:
        aeq_df.loc[:, cores] /= aeq_df[1]
    aeq_df.transpose().plot()
    aeq_df

cores | 1 | 2 | 4 | 8 | 16 | 32
--- | --- | --- | --- | --- | --- | ---
zones |  |  |  |  |  |
1500 | 1.0 | 0.559114 | 0.408111 | 0.250297 | 0.216308 | 0.255603
2500 | 1.0 | 0.526490 | 0.359411 | 0.255447 | 0.203943 | 0.232535
5000 | 1.0 | 0.517601 | 0.311059 | 0.216065 | 0.195720 | 0.187223
7500 | 1.0 | 0.528070 | 0.303304 | 0.221032 | 0.200886 | 0.189343
10000 | 1.0 | 0.520901 | 0.301453 | 0.209669 | 0.181450 | 0.174836
15000 | 1.0 | 0.520769 | 0.307402 | 0.208882 | 0.180139 | 0.176024