1
\$\begingroup\$

This question regards the C implementation of the python code posted in this other question of mine, hence the explanation of the task at hand will be the same:


The problem

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).


The question

I initially solved this problem with Python, but didn't manage to find a solution efficient enough. You can see the python code, and related discussions, on this codereview.SE question.

Here follows my solution implemented in C. This turned out to be much faster than the python solution. However, being me not an expert of C, I wonder if there is still room for improvements.

In particular, I'm wondering if it is better to read the whole file at once, which is what I've done in the code posted here, or to repeatedly read just a few bytes. In the linked question about the python implementation, switching from a read-all-at-once to a read-at-chunks solution gave a 2x improvement in speed. Doing the same thing with the C code, however, seems to significantly reduce the efficiency of the algorithm (though it can possibly be because of my poor way of implementing that).

Any tips/suggestions are welcome.


The code

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <sys/stat.h>
// #include <time.h>

off_t fsize(const char*);
int write_results_to_file(const char*, uint8_t(*)[4], uint16_t);

off_t fsize(const char* filename)
{
    struct stat st;

    if(stat(filename, &st) == 0)
        return st.st_size;

    return -1;
}

int write_results_to_file( const char* file_name,\
                           uint8_t (*results)[4], \
                           uint16_t number_of_elements
                          )
{
    FILE *file;
    file = fopen(file_name, "a");

    for(int i = 0; i < number_of_elements; i++)
    {
        printf("fourfold coincidence: ");
        for(int j = 0; j < 3; j++)
            printf("%d,", results[i][j]/32 + 1);
        printf("%d\n", results[i][3]/32 + 1);

        fprintf(file, "%d,%d,%d,%d\n",\
            results[i][0]/32 + 1,\
            results[i][1]/32 + 1,\
            results[i][2]/32 + 1,\
            results[i][3]/32 + 1
        );
    }

    fclose(file);
    return 0;
}

int main() {

    // clock_t begin_time;
    // begin_time = clock();

    const char FILE_NAME[] = "timetags.bin";
    const char OUTPUT_FILE_NAME[] = "processed_timetags.txt";
    const uint8_t WINDOW_SIZE = 100;
    // we use arrays with 4 elements, being interested only in fourfold coincidences. 
    uint32_t timestamps[4];
    uint8_t channels[4];
    // assuming that there are no more than 1000 coincidences
    uint8_t quadruples[1000][4];
    uint16_t quadruples_pos = 0;
    // curr_byte and curr_timestamp will be used in the for loop
    uint32_t curr_byte;
    uint32_t curr_timestamp;

    // Initialize array of chars able to contain all the bytes.
    // The first 40 bytes are to be ignored.
    // If the number of bytes is not a multiple of 5, the exceeding bytes are
    // not read.
    // In this way, file_size effectively equals 5 times the number of 5-byte blocks
    size_t file_size = fsize(FILE_NAME) - 40;
    file_size = file_size - file_size % 5;
    uint8_t* data = malloc(sizeof(uint8_t) * file_size);

    // open file in binary read mode
    FILE *file;
    file = fopen(FILE_NAME, "rb");
    if( file == NULL )
    {
        perror(FILE_NAME);
        return(-1);
    }
    // ignore the first 40 bytes
    fseek(file, 40, SEEK_SET);
    // the rest of the file is stored in data
    fread(data, file_size, 1, file);

    // We start by reading the first 5 bytes stored in data.
    // This is achieved reinterpreting the first 4 bytes in data as a single
    // 4-byte unsigned integer number.
    // The fifth byte encodes the channel.
    timestamps[0] = *(uint32_t*) &data[0];
    channels[0] = data[4];
    uint8_t timestamps_pos = 1;

    // Loop through data, at blocks of 5 bytes (remembering that we already
    // processed the first 5 bytes)
    for (int timetag_index = 1;
         timetag_index < file_size/5;
         timetag_index++
        )
    {
        // curr_byte points to the first byte of the 5-byte block that is
        // being processed
        curr_byte = timetag_index * 5;
        curr_timestamp = *(uint32_t*) &data[curr_byte];
        if( curr_timestamp - timestamps[0] <= WINDOW_SIZE &&
            curr_timestamp > timestamps[0]
          )
        {
            timestamps[timestamps_pos] = curr_timestamp;
            channels[timestamps_pos] = data[curr_byte + 4];
            timestamps_pos++;
        }
        else // we reinitialize the arrays and start a new window
        {
            // first the quadruple is stored, if we have one
            if(timestamps_pos == 4)
            {
                // quadruples[quadruples_pos++] = channels;
                // printf("found quadruple: ");
                // for(int i = 0; i<4; i++)
                    // printf("%d,\t",channels[i]);
                // printf("\t\tat timetag nr %d", timetag_index);
                // printf("\n");
                memcpy(&quadruples[quadruples_pos++], channels, sizeof(channels));
            }
            // then a new window is opened
            timestamps_pos = 1;
            timestamps[0] = curr_timestamp;
            channels[0] = data[curr_byte + 4];
        }
    }

    fclose(file);
    free(data);

    // printf("Time required: %f\n", (double) (clock()-begin_time) / CLOCKS_PER_SEC);

    write_results_to_file(OUTPUT_FILE_NAME, quadruples, quadruples_pos);

    return 0;

}
\$\endgroup\$
  • \$\begingroup\$ Reading from memory should always be faster unless the program itself (or other programs on your system) uses so much memory that you're using swap space. You might want to try the read-into-memory method again when your system is less loaded with processes. Obgrumble: back in MY day, we couldn't load 100M files into memory! \$\endgroup\$ – Barry Carter Jan 20 '16 at 15:07
  • 1
    \$\begingroup\$ For performance I would suggest to try memory-mapped-files. If your operating system does it's job properly you should have optimal performance while supporting files larger than memory. It does depend on the OS how to implement it, thus is not portable. \$\endgroup\$ – Zulan Jan 20 '16 at 17:15
  • \$\begingroup\$ @BarryCarter well, I'm running this program on a machine with 16GB of mamory, so that shuldn't be a problem. I really don't think it got to the point of having to use swap space. You can try for yourself to see if you get the same results: compare this version with this one \$\endgroup\$ – glS Jan 20 '16 at 19:18
1
\$\begingroup\$

Bug #1

You only leave room for 4 timestamps here:

uint32_t timestamps[4];

But if your file happens to contain 5+ timestamps in a row that all fit inside the time window, then your code will happily write past the end of the timestamps array. You need to check if you've reached 4 each time you add a new entry to the array.

When there are 5+ entries in a row that fit within the same window, I'm not sure if your program is supposed to output multiple quadruples or not, but you will need to figure out what you want to do in that case.

Bug #2

Your program currently can miss possible quadruples because it may skip past entries it should be looking at. For example, suppose the first six times were:

500, 599, 600, 601, 602, 900

with a time window of 100. Your program would start with time 500 and look for three more entries that fit within the 500..600 window. It would accept 599 and 600 but then fail at 601. But when it fails at 601, it will start a new window starting at 601 instead of "shifting the window by one" to start at 599. So it will end up missing the 599, 600, 601, 602 quadruple.

When you find an entry that falls outside the given window, you need to shift your window one by one until the newest entry falls inside the window. In other words, you throw out timestamps[0] and shift all the other entries down, and then retry the newest entry to see if it fits in the new window.

This could best be done with a circular buffer, where you wouldn't actually have to shift any entries (you would just shift the starting point forward).

Possible Bug #3

You mentioned that time could wrap around when it overflowed 4 bytes. But your code seems to reject any wrapped around times as not being valid:

    if( curr_timestamp - timestamps[0] <= WINDOW_SIZE &&
        curr_timestamp > timestamps[0]
      )

Perhaps that is what you intended, but if you wanted to find a quadruple that occurred near the wraparound time, you should allow wraparounds to be legal by removing the wraparound check like this:

    if (curr_timestamp - timestamps[0] <= WINDOW_SIZE)
\$\endgroup\$
  • \$\begingroup\$ thank you very much! I missed the first and third bug, but the second one is actually quite significant. I never heard of "circular buffer", I'll document myself and then see how to reimplement the code. \$\endgroup\$ – glS Jan 20 '16 at 19:26
  • \$\begingroup\$ @glS By "circular buffer", I mean you have an array of size 4, but instead of having array[0] always be the first entry, the first entry can be any of the 4 entries. So for example, after you read the first 4 entries, your starting point is [0] and the last entry is at [3]. But when you read the 5th entry, you read it into slot [0]. Now your starting point is [1] and your last entry is at [0] (wrapping around). \$\endgroup\$ – JS1 Jan 20 '16 at 19:39
3
\$\begingroup\$

Like @Barry Carter said in his comments, reading from memory should be much faster. Your algorithm looks pretty straightforward so I just have a few maintenance/style comments.

Remove Function Prototypes

In a single-file program like this, there's no need to have function prototypes sitting on top of the file. They only serve to give you another source of maintenance if your function signatures ever change.

Portability of fsize

Currently, your program will only work on a Linux computer. There's actually a standard library way of implementing fsize*:

size_t fsize(FILE *file)
{
    if (!file) { return 0; }

    const size_t curr_position = ftell(file);

    fseek(file, 0, SEEK_END);

    size_t file_size = ftell(file);

    fseek(file, curr_position, SEEK_SET);

    return size;
}

*NOTE: Although the standard mentions that SEEK_END is not required to be implemented, GCC, MSVC, and Clang support it.

Check for NULL after memory allocation

Your code never checks to see if the system was able to get the memory it asked for.

uint8_t* data = malloc(sizeof(uint8_t) * file_size);
if (!data) {
    fprintf(stderr, "Couldn't allocate %u bytes of data!\n", sizeof(uint8_t) * file_size);
    return 1;
}

Define a macro for 4

It seems like the number 4 is constantly used in your program as the size of arrays and such. I would put a macro definition for the number 4 on top of the source for maintainability.

#define TIMESTAMP_BYTE_COUNT 4

Miscellaneous

You input the wrong arguments to the fread function:

// the rest of the file is stored in data
    fread(data, file_size, 1, file);

fread takes a pointer to a buffer, the size of each element to be read in bytes, the number of elements to be read, and a file stream in that order. Also, don't put in a literal number when you can let the compiler figure out the size for you.

fread(data, sizeof(uint8_t), file_size, file);

I would separate this statement out by the way for readability:

memcpy(&quadruples[quadruples_pos++], channels, sizeof(channels));

into:

memcpy(&quadruples[quadruples_pos], channels, sizeof(channels));
++quadruples_pos;
\$\endgroup\$
  • \$\begingroup\$ thank you very much for pointing out these things. I actually run this program on a Windows machine, so that implementation of fsize does work on Windows, at least with MinGW's gcc. I see your point on portability though. Regarding the input arguments of fread, is there actually a difference in the way the two statements are implemented? the output seems to the be the same. \$\endgroup\$ – glS Jan 21 '16 at 8:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.