6
\$\begingroup\$

This is related to code I wrote for this SO question. Being most familiar with Python and Java, writing my first program in C was quite the learning curve. Any and all criticisms are welcome, especially anything I've handled unnaturally for the language.

Problem

Given a list of 1389 DNA strings of length 10 (complete file), e.g.,

ACTGCATGTC
CAACACAACG
TTCATGCCGA

I want to select 384 strings such that each letter in each position is as equally distributed as possible.

Approach

Represent each string with an array length 40, where each element describes whether a specific letter exists for each of the 10 positions. 1 for yes, 0 for no. Summing 384 such arrays should ideally result to 96 everywhere. Regardless of which strings are chosen, the mean of these summed arrays will always be 96, so we can computer the standard deviation (SD) to quantify how far away we're from an ideal solution.

We start by selecting any 384 random elements and keep the rest in a pool. For each of the chosen elements, we switch it with every one in the pool. If the SD improves, we make the switch. If no improvements are made after going through all 384 elements, the solution has converged. We can now feed these chosen indices back to the function, where it will switch one of the elements with one from the pool before converging this to a new solution. This feedback results to a biased sample of solutions with lower SD than if we were to choose new elements for each run independently.

Code

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>

#define NSTRINGS 1389
#define CHOOSE 384
#define LENGTH 40

void load_data(char *fname, int array[NSTRINGS][LENGTH]) {
    FILE *fp;
    char *line = NULL;
    size_t len = 0;
    ssize_t read;
    int counter = 0;

    fp = fopen(fname, "r");
    while ((read = getline(&line, &len, fp)) != -1) {
        for (int j = 0; j < strlen(line)-1; j++) {
            int stride;
            switch(line[j]) {
                case 'A':
                    stride = 0;
                    break;
                case 'T':
                    stride = 1;
                    break;
                case 'C':
                    stride = 2;
                    break;
                case 'G':
                    stride = 3;
                    break;
            }
            array[counter][10 * stride + j] = 1;
        }
        counter++;
    }
    fclose(fp);
}

void shuffle(int *array, size_t n) {
    for (int i = 0; i < n - 1; i++) {
        size_t j = i + rand() / (RAND_MAX / (n - i) + 1);
        int t = array[j];
        array[j] = array[i];
        array[i] = t;
    }
}

int element_in_array(int elem, int *array, size_t n) {
    for (int i = 0; i < n; i++) {
        if (elem == array[i]) {
            return 1;
        }
    }
    return 0;
}

void sum_rows(int table[NSTRINGS][LENGTH], int indices[CHOOSE], int sum[LENGTH]) {
    memset(sum, 0, sizeof(int) * LENGTH);
    for (int i = 0; i < CHOOSE; i++) {
        int *row = table[indices[i]];
        for (int j = 0; j < LENGTH; j++) {
            sum[j] += row[j];
        }
    }
}

double std(int array[LENGTH]) {
    double x = 0;
    double xsq = 0;

    for (int i = 0; i < LENGTH; i++) {
        int num = array[i];
        x += num;
        xsq += num * num;
    }

    xsq /= LENGTH;
    x /= LENGTH;

    return sqrt(xsq - x * x);
}

void init_arrays(int indices[CHOOSE], int pool[NSTRINGS-CHOOSE], int chain_solutions) {
    /*
     * Shuffling this array serves different purposes,
     * depending on which conditional block we enter.
     *
     * If it is to generate new indices, it's to choose
     * k elements from a population without replacement.
     *
     * If we already have an array of indices, it's to
     * choose a random index that isn't present already.
     */
    int idx[NSTRINGS];
    for (int i = 0; i < NSTRINGS; i++) {
        idx[i] = i;
    }
    shuffle(idx, NSTRINGS);

    if (!chain_solutions) {
        for (int i = 0; i < CHOOSE; i++) {
            indices[i] = idx[i];
        }
        for (int i = CHOOSE; i < NSTRINGS; i++) {
            pool[i-CHOOSE] = idx[i];
        }
    } else {
        /*
         * Algorithmically speaking, the code rotates out
         * the last element and considers the elements for
         * improvement in an ordered fashion. By shuffling
         * the indices array, we ensure there is no bias.
         */
        shuffle(indices, CHOOSE);
        for (int i = 0; i < NSTRINGS; i++) {
            if (!element_in_array(idx[i], indices, CHOOSE)) {
                indices[CHOOSE-1] = idx[i];
                break;
            }
        }
        for (int i = 0, j = 0; i < NSTRINGS; i++) {
            if (!element_in_array(i, indices, CHOOSE)) {
                pool[j++] = i;
            }
        }
    }
}

double minimise_variance(int table[NSTRINGS][LENGTH], int indices[CHOOSE], int chain_solutions) {
    int pool[NSTRINGS-CHOOSE];
    int temp_sum[LENGTH], final_sum[LENGTH];
    double metric;

    init_arrays(indices, pool, chain_solutions);
    sum_rows(table, indices, final_sum);
    metric = std(final_sum);
    while (1) {
        double start_metric = metric;
        for (int i = 0; i < CHOOSE; i++) {
            int index_row = indices[i];
            int pool_row = -1;
            /*
             * We want to collapse a (CHOOSE, LENGTH) array
             * to a (LENGTH) array by summing along the 0th
             * axis. Since CHOOSE-1 elements will always be
             * the same, we precompute their sum.
             */
            sum_rows(table, indices, temp_sum);
            for (int j = 0; j < LENGTH; j++) {
                temp_sum[j] -= table[index_row][j];
            }
            for (int j = 0; j < NSTRINGS-CHOOSE; j++) {
                for (int k = 0; k < LENGTH; k++) {
                    final_sum[k] = temp_sum[k] + table[pool[j]][k];
                }
                double measured_metric = std(final_sum);
                if (measured_metric < metric) {
                    metric = measured_metric;
                    pool_row = j;
                }
            }
            if (pool_row != -1) {
                indices[i] = pool[pool_row];
                pool[pool_row] = index_row;
            }
            else {
                indices[i] = index_row;
            }
        }
        if (start_metric == metric) {
            break;
        }
    }

    return metric;
}

int main(void) {
    // Commented out for the sake of repeatability for this review
    //srand(time(NULL));

    int runs = 10;
    int indices[CHOOSE] = {0};
    int table[NSTRINGS][LENGTH] = {0};
    load_data("final_list.txt", table);

    FILE *f = fopen("out", "wb");
    for (int i = 0; i < runs; i++) {
        printf("%d/%d\r", i+1, runs);
        minimise_variance(table, indices, i);
        fwrite(indices, sizeof(int), CHOOSE, f);
    }
    fclose(f);
}

Note

It wasn't worth the hassle, but if I wanted a more flexible program, I would have created a struct for the table and handled it like this.

typedef struct Table {
    size_t rows;
    size_t cols;
    int **arr;
} Table;

int main(void) {
    Table table;
    load_data(file, &table);
    // and later
    for (int i = 0; i < table.rows; i++) {
        free(table.arr[i]);
    }
    free(table.arr);
}

void load_data(char* fname, Table *table) {
    // read through the file once to get the number and length of strings, then
    int **arr = malloc(sizeof(*arr) * table->rows);
    for (int i = 0; i < table->rows; i++) {
        arr[i] = calloc(table->cols, sizeof(int));
    }
    table->arr = arr;
    // read data
}

Passing the various array lengths in functions would then be trivial without having to litter them with multiple size_t arguments. A pointer of pointers wouldn't have affected performance here because due to the nature of the code I never access sequential table rows anyway.

Regardless, I'm only interested in the actual code I've written, but if you think the hassle for a more dynamic approach is generally worth it, I'd like to hear it.

\$\endgroup\$
3
\$\begingroup\$

Very good first C code.

Sorry that most of the below is small stuff.

  1. load_data() should free(line) before returning. Note that getline() is not in the standard C library, yet a common extension.

     free(line); // add
     fclose(fp);
    
  2. Watch out when subtracting with unsigned math. strlen() returns type size_t, an unsigned type. With someone trying to hack code, the first character may be a null character, so strlen(line) == 0 and likely j < strlen(line)-1; will be the same as (size_t)j < SIZE_MAX, eventually causing int overflow with j++ and then undefined behavior.

    // Candidate solution, yet so below
    for (int j = 0; j + 1 < strlen(line); j++) {
    
  3. The above also implies various useful warnings are not enabled. Reccomend to enable all warnings.

  4. Assume the string length of line is N. A compiler may not recognize (or not allowed to given various reasons) that strlen(line) always givens the same answer as will instead call strlen(line) N times, - each time incurring a trip of O(N) as it finds the null character. So now the loop is O(n*n). Instead find the length and save it.

    size_t len = strlen(line);
    for (size j = 0; j + 1 < len; j++) {
    
  5. Yet even the above is not really needed. Just look for the end of the line.

    for (size_t j = 0; line[j] != '\n' && line[j] != '\0'; j++) {
    
  6. Use consistent types.

    // void shuffle(int *array, size_t n) {
    //    for (int i = 0; i < n - 1; i++) {        
    void shuffle(int *array, size_t n) {
        for (size_t i = 0; i + 1 < n; i++) {    
    
  7. Do not recommend the shuffle algorithm employed - it has weaknesses, Consider Fisher–Yates

  8. When possible, make array pointers point to const data. It self-documents code, allowed wider application and some optimizations.

    // int element_in_array(int elem, int *array, size_t n) {
    int element_in_array(int elem, const int *array, size_t n) {
    
    // double std(int array[LENGTH]) {
    double std(const int array[LENGTH]) {
    
  9. element_in_array() would make a good candidate to return bool rather than int.

  10. When code has the choice of the sizeof(some type) vs. sizeof some_object, choose the latter. Easier to code, review and maintain.

      // memset(sum, 0, sizeof(int) * LENGTH);
      memset(sum, 0, sizeof *sum * LENGTH);
      // or 
      memset(sum, 0, sizeof sum[0] * LENGTH);
    
    
      // fwrite(indices, sizeof(int), CHOOSE, f);
      fwrite(indices, sizeof *indices, CHOOSE, f);
    
  11. Beware of FP does not always compute the same as mathematical code would expect. xsq - x * x may be ever so slightly negative! Then sqrt(xsq - x * x); leads to unhappy code.

      // return sqrt(xsq - x * x);
      double sd2 = xsq - x * x;
      return sd2 < 0.0 ? 0.0 : sqrt(sd2); 
    
  12. Code is trying to minimize the standard deviation. Yet could instead minimize the square of the standard deviation. This is less computation (and slightly more accurate).

    double std2(const int array[LENGTH]) {
      double x = 0;
      ...
      // return sqrt(xsq - x * x);
      return xsq - x*x;
    }
    
  13. I'd expect more error checking on file I/O. Example:

    fp = fopen(fname, "r");
    // add
    if (fp == NULL) {
      Whine_and_complain(fname); 
      return;
    }
    
  14. A style issue, yet declaring and initializing is a nice coding style to consider. Avoids uninitialized variables and keeps them closer to use.

    // double metric;
    // metric = std(final_sum);
    double metric = std(final_sum);
    
\$\endgroup\$
  • \$\begingroup\$ This is exactly the feedback I was looking for! Subtle language idiosyncrasies that you may not come across the first time you're learning about the bigger picture. It's embarrassing I tripped up on 11 and 13. Regarding the bool type, I looked it up initially and saw it was effectively ints under the hood, which explains my decision. But I guess explicit use of bool makes my intentions clearer. \$\endgroup\$ – Reti43 Nov 15 '17 at 23:06
  • \$\begingroup\$ @Reti43 Since C99 _Bool/bool is its own type. It may or may not be the same space as int. Yet as you say, it does help convey code's intention. \$\endgroup\$ – chux Nov 15 '17 at 23:09
  • \$\begingroup\$ @Reti43 BTW: very good to have commented out //srand(time(NULL)); for review purposes. \$\endgroup\$ – chux Nov 15 '17 at 23:10
  • \$\begingroup\$ A clarification about #3. Are you referring to -Wall, or something else? Also, I have decided to change load_data() so that it returns -1 if either the file fails to open or an unrecognised character is encountered in the switch statement so that it behaves similarly to other reading functions. \$\endgroup\$ – Reti43 Nov 16 '17 at 1:05
  • 2
    \$\begingroup\$ @Reti43 I often use -Wall, -Wpedantic, -Wextra, -Wconversion with gcc. \$\endgroup\$ – chux Nov 16 '17 at 16:28
1
\$\begingroup\$

First it looks pretty good for a first time C program. Nice question! It is very nice to see a first C program without MAGIC Numbers.

The variable FILE *f in main() could probably use a more descriptive name.

The program could be more general. If the number of strings in the file increases the program will index past the end of the variable array. To make it more general, not waste space if there are less strings and not produce unknown results if there are more strings, the array could be implemented as a linked list. This would perform better than reading the file twice to get the size. The code never checks the value of counter to make sure it stays below 1388.

The code depends on fopen() working. There should be error checking for opening a file, reading from a file and writing to a file.

Technically there is nothing wrong with the switch statement in load_data(), but it would be better to have a default case that identifies what stride should if line[j] is not A, T, C or G.

        int stride;
        switch(line[j]) {
            case 'A':
                stride = 0;
                break;
            case 'T':
                stride = 1;
                break;
            case 'C':
                stride = 2;
                break;
            case 'G':
                stride = 3;
                break;
            default:
                stride = 0;
        }
\$\endgroup\$
  • 1
    \$\begingroup\$ Some good catches. I definitely fell into the mindset that things would work without failsafe checks because I was basing my implementation strictly to the file I had in mind. Not an excuse, but it is what it is. \$\endgroup\$ – Reti43 Nov 15 '17 at 23:09

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.