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.
