Efficiently processing large (~100 MB) structured binary data in Python 3

Question

I need to process the data contained in a (relatively) large binary file.

The file has the following structure:

40 bytes of initialization,
4 bytes, 1 byte,
4 bytes, 1 byte,
...
lots of such 5-byte blocks
...

The number of 5-byte blocks (to which I'll refer to as 'timetags' in the following) may vary, but the total size of the file can be in the order of ~100 MBs.

In each 5-byte block the first 4 bytes encode a uint32_t (unsigned integer) 'timestamp', and the fifth byte is a number encoding a 'channel'.

My task is to find out whether there are contiguous sequences of 4 timetags, such that the corresponding timestamps are within a certain time window from each other, and if that is the case store the corresponding channels.

For example, if there is a sequence of timetags whose decoded data is

100, 2
300, 4
310, 5
340, 8,
369, 6,
413, 8

and my time window is 100, then I store the list [4,5,8,6].

In general the number of such fourfold coincidences will be extremely small with respect to the total number of timetags (e.g. for a ~100MB file I have ~10 such coincidences). Also, the timestamps are generally in increasing order, but sometimes there is a sudden jump (when the timestamps becomes too big for the 4 bytes to encode) and the count starts over, and this has to be taken into account (see below for an example file).

I tried a variety of things to do this, but the best I achieved is a ~20s time to process a ~50MB file. However, implementing the same code in C, I get a performance of ~0.2s for a ~100MB file. While I understand that a C program is naturally faster than a python one, such a big performance gap really seems strange to me, so there must be something wrong with my python code.

Here is an example of the binary file to process (this is a reduced version of ~1MB). It contains a single fourfold coincidence at the last four timetags.

Here is my best implementation:

#!/usr/bin/python3
import struct
import os

# ============ CONSTANTS TO (OPTIONALLY) SET =============
WINDOW_SIZE = 100
INPUT_FILE_NAME = 'datafile.bin'

with open(INPUT_FILE_NAME, 'rb') as input_file:

    # number of bytes found in every file produced when writing the binary file
    initial_offset = 40
    # number of bytes occupied by each timetag
    timetag_size = 5

    # each timetag is stored in 5 bytes. The first 4 bytes contain the
    # timestamp, with the first byte being the least significant one, and last
    # byte the most significant one

    # bytes_to_timestamp expects in input a list of 4 bytes encoding a
    # timestamp. The first byte is the least significant one.
    def bytes_to_timestamp(bytes):
        return(struct.unpack("L", bytes)[0])

    # timestamps_in_window expects two 4-byte inputs, and returns whether the
    # second one should be considered in coincidence with the first one
    def timestamps_in_window(triggering_timestamp_bytes, new_timestamp_bytes):
        return(
                0 <
                bytes_to_timestamp(new_timestamp_bytes) -
                bytes_to_timestamp(triggering_timestamp_bytes)
                < WINDOW_SIZE
            )

    file_size = os.path.getsize(INPUT_FILE_NAME)
    effective_file_size = file_size - initial_offset - file_size % timetag_size

    input_file.seek(initial_offset)
    data = input_file.read(effective_file_size)

    output_channels = []
    triggering_timestamp = data[0:4]
    channels = [data[4]]
    for timetag_idx in range(effective_file_size // timetag_size - 1):
        current_byte = (timetag_idx + 1)*timetag_size
        current_timestamp = data[current_byte:current_byte+4]
        current_channel = data[current_byte + 4]
        if timestamps_in_window(triggering_timestamp, current_timestamp):
            channels.append(current_channel)
        else:
            if len(channels) == 4:
                print('found a fourfold coincidence at ', end = '')
                for channel in channels:
                    print('{0: >5}'.format(channel), end = '')
                print('')
                output_channels.append(channels)
            channels = [current_channel]
            triggering_timestamp = current_timestamp
    if len(channels) == 4:
        print('found a fourfold coincidence at ', end = '')
        for channel in channels:
            print('{0: >5}'.format(channel), end = '')
        print('')
        output_channels.append(channels)

I also tried to implement this with numpy.fromfile, but I ended up with a much slower code, which you can check out in this GitHub Gist.

How can I make the code more efficient?

Answer 1

Python isn't designed for raw speed. I'm not saying that it's slow, but that it makes tradeoffs in ways that don't always prioritize execution speed. Ease and expressiveness of coding tend to rate quite highly. So you have to learn some of the odd places that you can take better advantage of how Python is implemented.

Accessing variables; reducing indirection

When your code has a tight loop, C, for example, knows its variables' types, lays them out in memory and knows their addresses as it compiles, and can then directly manipulates bits. Python, largely to support its dynamic nature, has to do things a lot more indirectly. So it helps to know just how indirect some common things are. In particular, how indirect is looking up a name?

• If it's a local, CPython (usually) turns the name lookup into an array lookup
• If it's a global, it has to look up the name in the globals dictionary
• If it's an attribute or method, it has to look up the name in a sequence of dictionaries, depending on where it's first found, and, if relevant, a check and or call of a descriptor

When you're desperate to improve performance, try to get as many of these as possible to be the first case: looking up a local. In your original code, there is no function so nothing is a local. Thus step one is putting it in a function so most of your heavily used variables can be locals.

Sometimes you will want to consider further steps: giving a local name to a function or method that you call several times:

in_window = timestamps_in_window
append_channel = channels.append
: : :
for ...:
    if in_window(...):
        append_channel(...)

This avoids the overhead of looking up timestamps_in_window in globals(), or append in channels.__dict__ each iteration through the for loop, instead looking it up in the dictionary once, and ever after in the locals array.

Just be careful to avoid overdoing this. This sort of change can definitely harm readability (for example, what names would sufficiently help differentiate local names for channels.append and output_channels.append?), so this technique should be saved for critical performance paths, and its impact should be verified with representative profiling.

Reducing operations rather than strength

While strength reduction is useful in more direct languages like C because the processor is faster at certain kinds of operations, once you're in Python it's typically not as useful a metric. Instead you're better off focusing on reducing the number of operations, and knowing approximately each kind of operations relative expense. As a bonus, reducing the number of operations often makes code clearer.

So rather than trying to find tricks like x * 5 = (x << 2) + x, which may be faster in C, look for the opposite. Going the other way reduces two operations to one. And beware even the small x + 1 and x - 1 if you can remove them, as Python will have to call a method for each of these.

To make this concrete, consider these two lines of your code:

for timetag_idx in range(effective_file_size // timetag_size - 1):
    current_byte = (timetag_idx + 1)*timetag_size

You should be able to push the multiplication and division fully into the implementation layer by specifying a step for range:

for current_byte in range(0, effective_file_size, timetag_size):
    : : :

Then, to remove the + 1, replace the start (0) with timetag_size as well, or fold the initial case into your loop so you can start at 0.

That said, sometimes it does pay to think about the cost of a given operation. Unfortunately, like in many languages, function calls are relatively expensive. Unlike in many languages, there's little or no function inlining. This means there's a definite cost to clarifying names like bytes_to_timestamp. Why? This function means you're executing two function calls for a modest benefit in clarity.

Note that in vnp's answer, this function call is replaced with a single inlined call to struct.unpack that handles both the timestamp and channel in one call. It's possible that you could benefit from taking this further, and unrolling the loop a little bit more:

timetags = input_file.read(timetag.size * 4)
ts1, ch1, ts2, ch2, ts3, ch3, ts4, ch4 = struct.unpack("LBLBLBLB", timetags)

Be careful with this approach, however, as it will increase the complexity of the code up front, and then again once you have to deal with a file that includes timetags in quantities that aren't multiples of four. In general I would not expect this approach to be worth its cost, so definitely profile to prove whether it has value here if you use it.

Style

On the bright side, some of these speed recommendations line up with the style review of your code that I now offer. You have a good start; you comment your helper functions well. But you did not comment the longer blocks of (currently global) code that make up your program. Since this is typically the core of your algorithm, this may be a poor choice.

Putting them in functions will help in both ways. But as you put them in functions, be sure to consider what interface you should make available. Should the function know the filename, or should it be passed the open file object? Should it know what to do with coincident events, or should it be a generator that merely returns them as it finds them? I would tend towards the latter option in both cases:

def find_coincident_events(input_file, window_size):
    : : :
    if len(channels) == 4: # note others' comments about == vs >=
        yield channels

# global here, or could be a main()
with open(INPUT_FILE_NAME) as input_file:
    for coincidence in find_coincident_events(input_file, WINDOW_SIZE):
        print('found a fourfold coincidence at ', end = '')
        for channel in coincidence:
            : : :

Of course, even this separation might not be quite right. Perhaps the caller should seek past the header block, perhaps there are other parameters to consider. such as the a number of coincident events you require within a window, making it easy to search for 6 events in 120 units of time instead of 4 in 100.

Try to avoid duplication like you see in the two if len(channels) == 4: blocks. This could probably be refactored into a single function print_if_coincidence(channels), or at least the printing inside the if could be refactored into if len(channels) == 4: print_coincidence(channels) without seriously affecting performance (you say coincidences are rare). Note again the question of whether == or >= is correct for your problem statement, but realize that by putting it all into print_if_coincidence you would only have to update it in one place.

Answer 2

Reading the whole file at once is a major source of inefficiency. Think of all page faults and cache misses. Sometimes (for random access algorithms) it is unavoidable, but if the file is processed sequentially, streaming is more natural. You never need more than 2 records at a time, and you never seek backwards.

    try:
        while True:
            timetag = input_file.read(timetag.size)
            timestamp, channel = struct.unpack("LB", timetag)
            ....
    except struct.error:
        pass

Besides, your code recomputes each timestamp at least twice. The recommended approach does it just once, and simplifies comparison to

    def timestamps_in_window(triggering_timestamp_bytes, new_timestamp_bytes):
        return 0 < new_timestamp_bytes - triggering_timestamp_bytes < WINDOW_SIZE