7
\$\begingroup\$

This is an implementation of the viterbi algorithm in C, following from Durbin et. al.'s book Biological Sequence Analysis (2002). There's more info in the heading about usage and what exactle the program does.

It's working well and there are no memory leaks, but I'm interested in how I could structure the code better and perhaps make it run faster. I don't know much about good style in C yet. For one, I imagine there's too much in my main() function and I should break it down in to smaller functions, but I'm wondering the best way to do this.

I'm also interested in a better way to read in files, since I've heard that fscanf is risky and error prone. Any other error handling I'm missing as well?

If you want to compile and run the book example, there are the corresponding text files here: text files

/**
 * Implementation of the viterbi algorithm for estimating the states of a Hidden Markov Model given at least a sequence text file.
 * Program automatically determines n value from sequence file and assumes that state file has same n value.
 * 
 * Program follows example from Durbin et. al.'s "The occasionally dishonest casino, part 1." on p. 54 of Biological Sequence
 * Analyis (2002), with solution and viterbi output given on p. 57. The two states, F and L, correspond to a "Fair" or a "Loaded" die.
 * 
 * free .pdf of Durbin at: http://dnapunctuation.org/~poptsova/course/Durbin-Et-Al-Biological-Sequence-Analysis-CUP-2002-No-OCR.pdf
 * 
 * Optional argument to read in file of known states for comparison with algorithm's output. 
 * Sequence file and state files are is assumed to be one entry per line (see .txt files for example).
 * 
 * Usage: ./viterbi my_sequence_file.txt my_state_file.txt
 * my_sequence_file.txt = sequence file (required)
 * my_state_file.txt = state file (optional)
 *
**/


#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <ctype.h>

double max (double a, double b);
int argmax (double row0, double row1);

int main (int argc, char *argv[]) 
{
    // check for correct number of command line args
    if (argc != 2 && argc != 3)
    {
        printf("Usage: ./viterbi my_sequence_file.txt my_state_file.txt .  Include at least sequence file.\n");
        return 1;
    }

    // open sequence file and store in array. Dynamically allocate memory and automatically detect sequence n value
    FILE *seqf = fopen(argv[1], "r");
    if (!seqf)
    {
        printf("Invalid sequence file.\n");
        return 1;
    }

    int num;
    int memsize = 100;
    int n = 0;
    int *seq = calloc(memsize, sizeof(int));

    while(fscanf(seqf, "%i", &num) == 1)
    {
        seq[n] = num - 1;
        n++;

        if (n == memsize)
        {
            memsize += 100;
            seq = realloc(seq, memsize * sizeof(int));
        }

        if(!seq)
        {
            printf("Not enough memory.");
            return 1;
        }
    }
    fclose(seqf);

    // if passed as an argument, open the state solution file and print. Assumes n is same as sequence n above
    if(argv[2])
    {
        FILE *statef = fopen(argv[2], "r");

        if (!statef)
        {
            printf("Invalid state file.\n");
            return 1;
        }

        char *state = calloc(n, sizeof(char));
        if(!state)
        {
            printf("Not enough memory.");
            return 1;
        }

        char ch;
        printf("State solution:\n");
        for (int i = 0; i < n; i++)
        {
            fscanf(statef, "%c %*[\r\n]", &ch);
            state[i] = ch;

            if (i % 60 == 0 && i != 0)
            {
                printf("\n");
            }
            printf("%c", state[i]);
        }
        fclose(statef);
        free(state);
        printf("\n\n");
    }

    // state transition matrix in log space
    double a[2][2] = {
        { log(0.95),  log(0.05) },
        { log(0.1),  log(0.9) }
    };

    // emission probabilities, corresponding to p of rolling 1 thru 6 on fair or loaded die
    double e[6][2] = {
        { log( ((double) 1)/6),  log(0.1) },
        { log( ((double) 1)/6),  log(0.1) },
        { log( ((double) 1)/6),  log(0.1) },
        { log( ((double) 1)/6),  log(0.1) },
        { log( ((double) 1)/6),  log(0.1) },
        { log( ((double) 1)/6),  log(0.5) },
    };

    // allocate rest of memory and error handle
    int *path = calloc(n, sizeof(double));
    double **vprob = calloc(n, sizeof(double *));
    double **ptr = calloc(n, sizeof(double *));
    double **pi = calloc(n, sizeof(double *));

    if( !path || !vprob || !ptr || !pi )
        {
            printf("Not enough memory.");
            return 1;
        }

    for (int i = 0; i < 2; i++)
    {
        vprob[i] = calloc(n, sizeof(double));
        ptr[i] = calloc(n, sizeof(double));
        pi[i] = calloc(n, sizeof(double));

        if( !vprob[i] || !ptr[i] || !pi[i] )
        {
            printf("Not enough memory.");
            return 1;
        }
    }

    // initialize vprob array; assumed starting state is state F
    vprob[0][0] = 1;
    vprob[1][0] = 0;
    double row0;
    double row1;

    // viterbi algorithm in log space to avoid underflow
    for (int i = 1; i < n; i++)
    {
        for (int j = 0; j < 2; j++)
        {
            row0 = (vprob[0][i - 1] + a[0][j]);
            row1 = (vprob[1][i - 1] + a[1][j]);

            vprob[j][i] = e[seq[i]][j] + max( row0, row1 );
            ptr[j][i] = argmax( row0, row1 );
            pi[j][i] = max( row0 , row1 );
        }
    }
    free(seq);

    // traceback to find most likely path
    path[n - 1] = argmax( pi[0][n - 1], pi[1][n - 1] );
    for (int i = n - 2; i > 0; i--)
    {
        path[i] = ptr[path[i + 1]][i + 1];
    }

    // free remaining memory
    for (int i = 0; i < 2; i++)
    {
        free(vprob[i]);
        free(ptr[i]);
        free(pi[i]);
    }
    free(vprob);
    free(ptr);
    free(pi);

    // print viterbi result
    printf("Viterbi output:\n");
    for (int i = 0; i < n; i++)
    {   
        if (i % 60 == 0 && i != 0)
        {
            printf("\n");
        }

        if (path[i] == 0)
        {
            printf("F");
        }
        if (path[i] == 1)
        {
            printf("L");
        }
    }
    printf("\n");
    free(path);
    return 0;
}


double max (double a, double b)
{
    if (a > b)
    {
        return a;
    }
    else if (a < b)
    {
    return b;
    }
    // if equal, returns arbitrary argument for specific use in this algorithm
    return b;
}

int argmax (double row0, double row1)
{
    if (row0 > row1)
    {
        return 0;
    }
    else if (row0 < row1)
    {
        return 1;
    }
    return row1;
}
\$\endgroup\$

4 Answers 4

6
\$\begingroup\$

Well, you are right that you rediscovered the megamoth, and it should be divided into easier to understand and maybe reusable parts.

But let's first look at your 2 functions, which are both unneccessarily long and verbose.

  1. If you just used else instead of else if you could subsume the third case (equality) with one of the other two.
  2. There is such a thing called the conditional operator, C's only ternary operator, for selecting one of two values.
  3. It looks like argmax() should return 1 in the equal case, instead of converting the first floating-point argument to integer, parallel to max().

That gives you:

double max(double a, double b)
{
    return a > b ? a : b;
}

int argmax(double a, double b)
{
    return a > b ? 0 : 1;
}
  1. You really should restrict yourself to some reasonable line-length like 79, horizontal scrolling is absolute murder.
    To that end, know that two string-literals put next to each other will be fused by the compiler:

    strcmp("a" "b", "ab") == 0
    
  2. Many of your comments restate the obvious. Others can be made superfluous by extracting a properly named function.
    Candidates are read_sequencefile(), read_and_echo_statefile(), do_analysis and print_result.

  3. There's a standard way to denote optional arguments in a usage-message: square brackets.
  4. There's no need to hardcode your executables filename, look at argv[0].
  5. If you use calloc you should actually benefit from getting the returned memory zeroed. Otherwise it's just useless makework.
  6. Try to minimize the scope of your variables. While the compiler can be expected to defer reservation and initialization and bring forward oblivion, that doesn't help you. The less there is to keep track of, the better.
  7. seq = realloc(seq, memsize * sizeof(int)); is not an error because you immediately exit the program. Otherwise, it would be a memory-leak.
  8. Allocate memory when you need it, not if you might need it later.
  9. Errors should go to stderr, and all output should end with a newline.
  10. If you only want to output a single char, consider using putc.
  11. You can use sizeof with expressions and types. Try to avoid the latter as it forces you to repeat your types.
  12. You should check whether your State file is corrupted.
  13. Anyway, did you ever want to do something with the State file?
  14. If you want to select something with an index and your indices are dense, consider using an array:

    putc("FL"[path[i]], stdout);
    // Equivalent if path[i] in 0, 1
    if (path[i] == 0)
    {
        printf("F");
    }
    if (path[i] == 1)
    {
        printf("L");
    }
    
  15. Instead of casting an integer-constant to double, consider simply making it a double-constant: 1.
  16. If you allocate and deallocate multiple pieces at the same times, consider coalescing them to only bother the allocator once.
  17. Only the first two elements of vprob, ptr and pi are ever used. Don't allocate more, and allocate such small arrays on the stack.
  18. Mark anything which shall not be modified const. That frees everyone of the fear of missing a change, and lets the compiler enforce it.
  19. return 0; is implied for main since C99.
  20. EXIT_SUCCESS resp. 0 and EXIT_FAILURE are the only standardized exit-codes.
  21. Consider making all functions only used in the local file local with static. Doing so encourages inlining and avoids emitting a symbol for the linker.

Refactored program:

/**
 * Implementation of the viterbi algorithm for estimating the states of a
 * Hidden Markov Model given at least a sequence text file.
 * Program automatically determines n value from sequence file and assumes that
 * state file has same n value.
 * 
 * Program follows example from Durbin et. al.'s "The occasionally dishonest
 * casino, part 1." on p. 54 of Biological Sequence Analyis (2002), with
 * solution and viterbi output given on p. 57. The two states, F and L,
 * correspond to a "Fair" or a "Loaded" die.
 * 
 * free .pdf of Durbin at:
 * http://dnapunctuation.org/~poptsova/course/Durbin-Et-Al-Biological-Sequence-Analysis-CUP-2002-No-OCR.pdf
 * 
 * Optional argument to read in file of known states for comparison with
 * algorithm's output.
 * Sequence file and state files are is assumed to be one entry per line (see
 * .txt files for example).
 * 
 * Usage: ./viterbi my_sequence_file.txt my_state_file.txt
 * my_sequence_file.txt = sequence file (required)
 * my_state_file.txt = state file (optional)
 *
**/


#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <ctype.h>

static double max(double a, double b)
{
    return a > b ? a : b;
}

static int argmax(double a, double b)
{
    return a > b ? 0 : 1;
}

static int read_sequencefile(const char filename[], int** sequence)
{
    assert(sequence);
    FILE* file = fopen(filename, "rb");
    if (!file) {
        fprintf(stderr, "Invalid sequence file.\n");
        exit(EXIT_FAILURE);
    }

    int size = 0;
    int* seq = 0;

    for(int num, capacity = 0; fscanf(file, "%i", &num) == 1; ) {
        if(size == capacity) {
            capacity += 100;
            seq = realloc(seq, capacity * sizeof *seq);
            if(!seq) {
                fprintf(stderr, "Not enough memory.\n");
                exit(EXIT_FAILURE);
            }
        }
        seq[size++] = num;
    }
    if(!feof(file)) {
        fprintf(stderr, "Sequence file contains garbage.\n");
        exit(EXIT_FAILURE);
    }
    fclose(file);
    *sequence = seq;
    return size;
}

static void read_and_echo_statefile(const char filename[], int n)
{
    FILE *file = fopen(argv[2], "r");
    if (!file) {
        fprintf(stderr, "Invalid state file.\n");
        exit(EXIT_FAILURE);
    }

    char* state = malloc(n * sizeof *state);
    if(n && !state) {
        fprintf(stderr, "Not enough memory.\n");
        exit(EXIT_FAILURE);
    }

    printf("State solution:\n");
    for (int i = 0; i < n; i++) {
        char ch;
        if(fscanf(file, " %c ", &ch) != 1) {
            fprintf(stderr, "State file too small.\n");
            exit(EXIT_FAILURE);
        }
        state[i] = ch;

        if (i % 60 == 0 && i != 0)
            putc('\n', stdout);
        putc(state[i], stdout);
    }
    if(!feof(file)) {
        fprintf(stderr, "State file too large.\n");
        exit(EXIT_FAILURE);
    }
    free(state);
    fclose(file);
    printf("\n\n");
}

static int* do_analysis(const int seq[], int n)
{
    // state transition matrix in log space
    const double a[2][2] = {
        { log(0.95),  log(0.05) },
        { log(0.1),  log(0.9) }
    };

    // emission probabilities, corresponding to probability of rolling 1
    // through 6 on fair or loaded die
    const double e[6][2] = {
        { log( 1. / 6),  log(0.1) },
        { log( 1. / 6),  log(0.1) },
        { log( 1. / 6),  log(0.1) },
        { log( 1. / 6),  log(0.1) },
        { log( 1. / 6),  log(0.1) },
        { log( 1. / 6),  log(0.5) },
    };

    int *path = malloc(n * sizeof *path);
    *path = 0;
    double* vprob[2] = { calloc(3 * 2 * n, sizeof **vprob), NULL };
    if(n && (!path || ! vprob[0])) {
        fprintf(stderr, "Not enough memory.\n");
        exit(EXIT_FAILURE);
    }
    vprob[1] = vprob[0] + n;
    double* ptr[2] = { vprob[1] + n, vprob[1] + 2 * n };
    double* pi[2] = { ptr[1] + n, ptr[1] + 2 * n };

    // initialize vprob array; assumed starting state is state F
    vprob[0][0] = 1;
    vprob[1][0] = 0;

    // viterbi algorithm in log space to avoid underflow
    for (int i = 1; i < n; i++)
    {
        for (int j = 0; j < 2; j++)
        {
            double row0 = (vprob[0][i - 1] + a[0][j]);
            double row1 = (vprob[1][i - 1] + a[1][j]);

            vprob[j][i] = e[seq[i]][j] + max( row0, row1 );
            ptr[j][i] = argmax( row0, row1 );
            pi[j][i] = max( row0 , row1 );
        }
    }

    // traceback to find most likely path
    path[n - 1] = argmax( pi[0][n - 1], pi[1][n - 1] );
    for (int i = n - 2; i > 0; i--)
        path[i] = ptr[path[i + 1]][i + 1];

    free(vprob[0]);
    return path;
}

static void print_result(const int path[], int n)
{
    printf("Viterbi output:\n");
    for (int i = 0; i < n; i++) {   
        if (i % 60 == 0 && i != 0)
            putc('\n', stdout);
        putc("FL"[path[i]], stdout);
    }
    putc('\n', stdout);
}

int main(int argc, char *argv[]) 
{
    if (argc != 2 && argc != 3) {
        printf("Usage: %s sequence_file.txt [state_file.txt]\n", argv[0]);
        return EXIT_FAILURE;
    }

    int* seq;
    int n = read_sequencefile(argv[1], &seq);

    if(argv[2])
        read_and_echo_statefile(argv[2], n);

    int* path = do_analysis(seq, n);
    free(seq);

    print_result(path, n);
    free(path);
}

I didn't test all occurrences of calloc whether malloc is sufficient.
The only line 80 or bigger contains the URL, because you just don't break URLs if possible.

\$\endgroup\$
4
  • \$\begingroup\$ A lot of great stuff to sift through here between your post and the one below. Thanks so much. \$\endgroup\$
    – Davigor
    May 18, 2017 at 22:09
  • \$\begingroup\$ Very complete review! Just a very small spot about another latent error: you assign to int* the memory allocated for sizeof(double). I'm not aware of any system where sizeof(int) > sizeof(double) but in that case it will fail when you deference that pointer. \$\endgroup\$ May 19, 2017 at 8:03
  • \$\begingroup\$ Good deep review. 1) double* vprob[2] = { calloc(3 * 2 * n, sizeof **vprob) }; looks wrong with only 1 initializer. Suggest should be double* vprob[2] = { calloc(3 * 2 * n, sizeof *vprob[0]), NULL };` 2) An empty file will int n = 0; and could take the if path with if(!path || ! vprob[0]) { provide erroneous "Not enough memory.\n" message. 3) Per #21, void print_result(int* path, int n) --> void print_result(const int* path, int n) 4) putc("FL"[path[i]], stdout); treads on thin ice. Perhaps putc("FL"[1&path[i]], stdout); or [(bool)path[i]] or use bool* path \$\endgroup\$ May 19, 2017 at 16:45
  • \$\begingroup\$ ... 5) C does not specified "integer-literal" but integer constant. Like-wise for `double literal. The 2 literals C specifies (string, compound) can have their address taken unlike an integer constant. \$\endgroup\$ May 19, 2017 at 16:48
3
\$\begingroup\$

Welcome to Code Review, very nice first question for Code Review!

You're quite correct about the program structure and you may want to do some research on software design principles such as the Single Responsibility Principle, the Don't Repeat Yourself Principle (DRY code), the Keep It Simple KIS(S) principle and Demeter's Law.

What is the State File For?

Currently the program is reading the state information from the file, but the state information is never used, the state array is thrown away after it is input. This could almost be considered a bug or broken code.

Use stderr to Report Errors

When C was first developed the creators created 3 File pointers to use, stdin, stdout and stderr. The function printf() prints to stdout, scanf() reads from stdin and a special file pointer for reporting error messages was created, this is stderr. It is better to report errors to stderr than stdout because redirecting file output by progname arguments > outfile does not redirect stderr. It is possible to redirect stderr to a file as well using progname arguments >& outfile (Unix and Linux).

To print errors to stderr:

    fprintf(stderr, "ERROR MESSAGE");

Reduce Complexity, Follow SRP

The Single Responsibility Principle states that every module or class should have responsibility over a single part of the functionality provided by the software, and that responsibility should be entirely encapsulated by the class. All its services should be narrowly aligned with that responsibility.

Robert C. Martin expresses the principle as follows:
    `A class should have only one reason to change.`

While this is primarily targeted at classes in object oriented languages it applies to functions and subroutines in procedural languages like C as well.

The main function could be broken up into at least 4 functions:

int* resize_seq(int* seq, int* memsize);
int* process_sequence_input_file(/* in variable */ char* input_file, /* out variables */ int* out_n)
char* process_state_file(char* statefilename, int n);
int execute_biological_sequence_analysis(int *seq, int sequence_size, char *state)

The function execute_biological_sequence_analysis() can be broken up into multiple functions as well.

Here is an example of what the program might look like following the Single Responsibility Principle.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <ctype.h>

#define SEQUENCE_MEMORY_ALLOCATION_SIZE    100
#define DIE_FACE_COUNT  6

double max (double a, double b)
{
    if (a > b)
    {
        return a;
    }
    else if (a < b)
    {
        return b;
    }
    // if equal, returns arbitrary argument for specific use in this algorithm
    return b;
}

double argmax (double row0, double row1)
{
    if (row0 > row1)
    {
        return 0;
    }
    else if (row0 < row1)
    {
        return 1;
    }
    return row1;
}

int* resize_seq(int* seq, int* memsize)
{
    *memsize += SEQUENCE_MEMORY_ALLOCATION_SIZE;
    seq = realloc(seq, *memsize * sizeof(int));
    if(!seq)
    {
        fprintf(stderr, "Not enough memory to allocate or reallocate seq.");
    }
    return seq;
}

int* process_sequence_input_file(/* in variable */ char* input_file, /* out variables */ int* out_n)
{
    int memsize = SEQUENCE_MEMORY_ALLOCATION_SIZE;
    int n = 0;
    int num;

    int *seq = calloc(memsize, sizeof(int));
    if (!seq)
    {
        return NULL;
    }

    // open sequence file and store in array. Dynamically allocate memory and automatically detect sequence n value
    FILE *seqf = fopen(input_file, "r");
    if (!seqf)
    {
        fprintf(stderr, "Can't open input file %s.\n", input_file);
        free(seq);
        return NULL;
    }

    while(fscanf(seqf, "%i", &num) == 1)
    {
        seq[n] = num - 1;
        n++;

        if (n == memsize)
        {
            if (!(seq = resize_seq(seq, &memsize))) {
                return NULL;
            }
        }

    }
    fclose(seqf);

    *out_n = n;         // return n

    return seq;
}

char* process_state_file(char* statefilename, int n)
{
    FILE *statef = fopen(statefilename, "r");

    if (!statef)
    {
        printf("Invalid state file.\n");
        return NULL;
    }

    char *state = calloc(n, sizeof(char));
    if(!state)
    {
        printf("Not enough memory.");
        return NULL;
    }

    char ch;
    printf("State solution:\n");
    for (int i = 0; i < n; i++)
    {
        fscanf(statef, "%c %*[\r\n]", &ch);
        state[i] = ch;

        if (i % 60 == 0 && i != 0)
        {
            printf("\n");
        }
        printf("%c", state[i]);
    }
    fclose(statef);
    printf("\n\n");

    return state;
}

int execute_biological_sequence_analysis(int *seq, int sequence_size, char *state)
{
    int status = EXIT_SUCCESS;

    // state transition matrix in log space
    double a[2][2] = {
        { log(0.95),  log(0.05) },
        { log(0.1),  log(0.9) }
    };

    // emission probabilities, corresponding to p of rolling 1 thru 6 on fair or loaded die
    double e[DIE_FACE_COUNT][2] = {
        { log( ((double) 1)/DIE_FACE_COUNT),  log(0.1) },
        { log( ((double) 1)/DIE_FACE_COUNT),  log(0.1) },
        { log( ((double) 1)/DIE_FACE_COUNT),  log(0.1) },
        { log( ((double) 1)/DIE_FACE_COUNT),  log(0.1) },
        { log( ((double) 1)/DIE_FACE_COUNT),  log(0.1) },
        { log( ((double) 1)/DIE_FACE_COUNT),  log(0.5) },
    };

    // allocate rest of memory and error handle
    int *path = calloc(sequence_size, sizeof(double));
    double **vprob = calloc(sequence_size, sizeof(double *));
    double **ptr = calloc(sequence_size, sizeof(double *));
    double **pi = calloc(sequence_size, sizeof(double *));

    if( !path || !vprob || !ptr || !pi )
    {
        printf("Not enough memory.");
        return 1;
    }

    for (int i = 0; i < 2; i++)
    {
        vprob[i] = calloc(sequence_size, sizeof(double));
        ptr[i] = calloc(sequence_size, sizeof(double));
        pi[i] = calloc(sequence_size, sizeof(double));

        if( !vprob[i] || !ptr[i] || !pi[i] )
        {
            printf("Not enough memory.");
            return EXIT_FAILURE;
        }
    }

    // initialize vprob array; assumed starting state is state F
    vprob[0][0] = 1;
    vprob[1][0] = 0;
    double row0;
    double row1;

    // viterbi algorithm in log space to avoid underflow
    for (int i = 1; i < sequence_size; i++)
    {
        for (int j = 0; j < 2; j++)
        {
            row0 = (vprob[0][i - 1] + a[0][j]);
            row1 = (vprob[1][i - 1] + a[1][j]);

            vprob[j][i] = e[seq[i]][j] + max( row0, row1 );
            ptr[j][i] = argmax( row0, row1 );
            pi[j][i] = max( row0 , row1 );
        }
    }

    // traceback to find most likely path
    path[sequence_size - 1] = argmax( pi[0][sequence_size - 1], pi[1][sequence_size - 1] );
    for (int i = sequence_size - 2; i > 0; i--)
    {
        path[i] = ptr[path[i + 1]][i + 1];
    }

    // free remaining memory
    for (int i = 0; i < 2; i++)
    {
        free(vprob[i]);
        free(ptr[i]);
        free(pi[i]);
    }
    free(vprob);
    free(ptr);
    free(pi);

    // print viterbi result
    printf("Viterbi output:\n");
    for (int i = 0; i < sequence_size; i++)
    {
        if (i % 60 == 0 && i != 0)
        {
            printf("\n");
        }

        if (path[i] == 0)
        {
            printf("F");
        }
        if (path[i] == 1)
        {
            printf("L");
        }
    }
    printf("\n");
    free(path);

    return status;
}

int main (int argc, char *argv[])
{
    int status = EXIT_SUCCESS;
    // check for correct number of command line args
    if (argc != 2 && argc != 3)
    {
        fprintf(stderr, "Usage: ./viterbi my_sequence_file.txt my_state_file.txt .  Include at least sequence file.\n");
        return EXIT_FAILURE;
    }

    int sequence_size = 0;
    int *seq;

    if (!(seq = process_sequence_input_file(argv[1], &sequence_size)))
    {
        return EXIT_FAILURE;
    }

    char *state = NULL;
    // if passed as an argument, open the state solution file and print. Assumes n is same as sequence n above
    if(argv[2])
    {
        if (!(state = process_state_file(argv[2], sequence_size)))
        {
            return EXIT_FAILURE;
        }
    }
    status = execute_biological_sequence_analysis(seq, sequence_size, state);

    free(seq);
    free(state);

    return status;
}

Use Sybolic Constants Rather than Numeric Constants

There are a number of numeric constants such as 6, 100, 0.1, 0.5 that could be represented as symbolic constants, here is an example of two of them

#define SEQUENCE_MEMORY_ALLOCATION_SIZE    100
#define DIE_FACE_COUNT  6

There are 2 constants defined in stdlib.h and stdlib, these are EXIT_SUCCESS and EXIT_FAILURE. These make main() more readable rather than return 0, return EXIT_SUCCESS or EXIT_FAILURE from main, this will also allow error handling to be added at some point in the future.

Symbolic constants help programmers in a number of ways:

  • They make the code more readable
  • They decrease the number of lines that need to be changed in an edit, this makes changing the code easier in the future.
  • When using constants to create the size of arrays, it is much easier to change all the arrays sizes and loop constants at once.

Variable Names

Generally the variable names are ok, some might improve with additional length such as changing seq to sequence. One variable should probably be renamed from n to `sequence_size since it is used in much of the program.

Testing the Return of Allocation

The code contains all the necessary testing to see if calloc() returns a value or not, but it would be better if each call to calloc() was immediately tested rather than what the following code does:

    // allocate rest of memory and error handle
    int *path = calloc(sequence_size, sizeof(double));
    double **vprob = calloc(sequence_size, sizeof(double *));
    double **ptr = calloc(sequence_size, sizeof(double *));
    double **pi = calloc(sequence_size, sizeof(double *));

    if( !path || !vprob || !ptr || !pi )
    {
        printf("Not enough memory.");
        return EXIT_FAILURE;
    }

    for (int i = 0; i < 2; i++)
    {
        vprob[i] = calloc(sequence_size, sizeof(double));
        ptr[i] = calloc(sequence_size, sizeof(double));
        pi[i] = calloc(sequence_size, sizeof(double));

        if( !vprob[i] || !ptr[i] || !pi[i] )
        {
            printf("Not enough memory.");
            return EXIT_FAILURE;
        }
    }

One possible solution would be to create a function that calls calloc() and then tests the return value, this would apply the DRY principle (Don't Repeat Yourself).

double** calloc_with_test_and_error_message(int memsize, char *arrayName)
{
    double** arrayofdoublepointers = calloc(1, memsize);
    if (!arrayofdoublepointers)
    {
         fprintf(stderr, "ERROR: calloc failed to create the array of double pointers %s\n", arrayName);
    }
    return arrayofdoublepointers;
}
\$\endgroup\$
7
  • \$\begingroup\$ Thanks for such a comprehensive reply, a lot of goof stuff for me to check out here. \$\endgroup\$
    – Davigor
    May 18, 2017 at 22:09
  • \$\begingroup\$ int argmax (double row0, double row1) { ... return row1; } looks wrong. Returning a double as an int? \$\endgroup\$ May 19, 2017 at 16:50
  • \$\begingroup\$ int* process_sequence_input_file(/* in variable */ char* input_file, ... is sufficient as int* process_sequence_input_file(const char* input_file, ... as the const enforces the "in-ness" of input_file. ` \$\endgroup\$ May 19, 2017 at 16:52
  • \$\begingroup\$ if(!state) { printf("Not enough memory."); is not necessarily true if n==0. Could use if(!state && n).... IMO, better for n to be size_t than int, yet that reflects back on much of OP code. \$\endgroup\$ May 19, 2017 at 16:56
  • \$\begingroup\$ Using ptr = alloc(sizeof *ptr * n) model would have not caused the coding oops in int *path = calloc(sequence_size, sizeof(double));. Recommend int *path = calloc(sequence_size, sizeof *path); \$\endgroup\$ May 19, 2017 at 17:00
2
\$\begingroup\$

Too large allocations

I noticed that you allocated several arrays of length n here:

// allocate rest of memory and error handle
int *path = calloc(n, sizeof(double));
double **vprob = calloc(n, sizeof(double *));
double **ptr = calloc(n, sizeof(double *));
double **pi = calloc(n, sizeof(double *));

if( !path || !vprob || !ptr || !pi )
    {
        printf("Not enough memory.");
        return 1;
    }

Problem number one is that for path you allocated space for n doubles instead of n ints (typo most likely). To avoid this problem you should use the variable itself in the sizeof() instead of hardcoding a type.

Problem number two is that of these 4 arrays, only path needed to be of size n. The other three are only supposed to be size 2. So really, you don't even need to allocate those, you could just do this:

// allocate rest of memory and error handle
int    *path    = calloc(n, sizeof(*path));
double *vrob[2] = {0};
double *ptr [2] = {0};
double *pi  [2] = {0};

if( !path )
{
    printf("Not enough memory.");
    return EXIT_FAILURE;
}

Potential buffer overflow

Your argmax() function returns row1 when row0 == row1, but I think it is supposed to return 1 instead. Later on, the result of argmax() is used as an index into the ptr[] array, which only has two entries. So if row0 and row1 were both 5, for example, you would get a buffer overflow later on.

\$\endgroup\$
1
  • \$\begingroup\$ Idea: As vrob[2] is fixed at 2, double *vrob[2] = {NULL, NULL}; adds clarity. if( !path )has same corner case issue when n==0 in providing wrong message. \$\endgroup\$ May 19, 2017 at 17:05
2
\$\begingroup\$

OP's code is largely well reviewed. Just adding a bit


The "%*[\r\n]" in fscanf(statef, "%c %*[\r\n]", &ch); serves no purpose.

The " " in fscanf() format consumes any 0 or more white-space characters including \r and \n. So these is nothing left to consume by "%*[\r\n]".

Further, if the file begins with a white space like '\n', that character will get saved in ch. Not likely code's intent.

If there is not enough data, code could be reading uninitialized data with state[i] = ch;, which is undefined behavior (UB).

Consider validating fscanf() results before using ch by checking, at least, the return value.

    char ch;
    printf("State solution:\n");
    fflush(stdout);  // Insure output is seen 
    for (int i = 0; i < n; i++) {
        // skip white-spaces 
        //                  v
        if (fscanf(statef, " %c", &ch) != 1) {
          Handle_Scant_Input();
        }
        state[i] = ch;
        ...
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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