evsim.c version 0.1

// evsim.c
// Simple simulation of genetic drift
// Arguments: m n g l s d
// m = mutation rate (capped at 1.0) per birth
// n = population size [default: 1000]
// g = generations (needs to be at least 4*n to get close to steady state) [default: 40*n]
// l = fraction of mutatinos that are lethal (0.01 = 1%) [default: 0.01]
// s = fraction of mutations that are strongly negative (10% lower survival chance) [default: 0.02]
// d = fraction of mutations that are mildly negative (1% lower survival chance) [default: 0.03]
//
// If desired, consider a 'p' that is fraction of mutations that increase survival by 1%
//
// Version 0.1 - David Fenger

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

#define GENE_BYTES 10000

float m_rate;
long pop_size = 1000;
long generations, cur_generation;
long cur_pop;
long min_pop;
float neg_l=0.01,neg_s=0.02,neg_m=0.03;
long neg_lb,neg_sb,neg_mb;
unsigned char *p_genes, *p_offspring;
unsigned char fixed_genes[GENE_BYTES];
long c_fixed = 0;
long first_fix = 0;
long first_half_fixed = -1;
float inv_rmax = 1.0f / (float)RAND_MAX;
char bit_count[256];

// TODO: Find a good rand(), don't forget to seed it.
// TODO: Print seed value used, for future reference.  And add as another argument...  [default = clock]

float frand() {
    return (float)rand() * inv_rmax;
}

long lrand(long n) {
    return rand() % n;
}

// Note: this is not as random as I'd like, as even a 1% mutation rate will occasionally
// cause a double mutation (1 time in 10,000).  But it's close enough for the purpose, and fast.
void copy_with_mutate(int icur, int ioff) {
    long i,loc;

    memcpy(p_offspring+ioff*GENE_BYTES, p_genes+icur*GENE_BYTES, GENE_BYTES);

    // Roll for mutation
    if(frand() <= m_rate) {
        loc = lrand(GENE_BYTES*8);
        p_offspring[ioff*GENE_BYTES+loc/8] |= 1<<(loc%8);
    }
}

void copy_to_pop(long icur, long ioff)
{
    memcpy(p_genes+icur*GENE_BYTES, p_offspring+ioff*GENE_BYTES, GENE_BYTES);
}

void do_breed() {
    // Fill offspring array
    long i;
    // To be fair, everyone gets at least one offspring
    for(i=0;i<cur_pop;i++) {
        copy_with_mutate(i,i);
    }
    // Need 2xpop_size offspring to cull from
    for(;i<2*pop_size;i++) {
        copy_with_mutate(lrand(cur_pop), i);
    }
}

float calc_viability(long ioff)
{
    float v = 0.5f; // Baseline survivability
    long i;
    char *p = p_offspring+ioff*GENE_BYTES;
    // Lethality - return 0 for any lethal mutations
    for(i=0;i<neg_lb;i++) {
        if(p[i] > 0) return 0.0f;
    }

    // Strongly negative mutations - each reduces v by 0.1;
    for(;i<neg_lb+neg_sb;i++) {
        if(p[i] > 0) v -= bit_count[p[i]] * 0.1;
    }
    if(v <= 0.0f) return 0.0f;

    // Mildly negative mutations - each reduces v by 0.01;
    for(;i<neg_lb+neg_sb+neg_mb;i++) {
        if(p[i] > 0) v -= bit_count[p[i]] * 0.01;
    }
    if(v <= 0.0f) return 0.0f;
    return v;
}

void do_cull() {
    // Kill off children, copying survivors into original population
    long ipop = 0;
    long ioff = 0;
    float v,vt;
    vt = 0;
    for(ioff = 0; ioff < 2*pop_size; ioff++) {
        // Compute viability of offspring
        v = calc_viability(ioff);
        vt+=v;
        if(frand() < v) {
            copy_to_pop(ipop, ioff);
            ipop++;
        }
    }
    cur_pop = ipop;
    if(cur_pop < min_pop) min_pop = cur_pop;
    if(cur_pop <= 0) {
        printf("WHOOPS - population is dead.\n");
        return(-5);
    }
}

void do_evaluate(long cur_generation) {
    // Count # of fixed genes
    long ipop, igene, c;

    memcpy(fixed_genes, p_genes, GENE_BYTES); // Seed with first organism

    for(ipop = 1; ipop < cur_pop; ipop++) {
        for(igene=0;igene<GENE_BYTES;igene++) {
            fixed_genes[igene] &= p_genes[ipop*GENE_BYTES+igene];
        }
    }

    // Count bits set in fixed_genes - these are mutations possessed by entire population
    c = 0;
    for(igene=0;igene<GENE_BYTES;igene++) {
        c += bit_count[fixed_genes[igene]];
    }
    if(c > c_fixed) {
        if(c_fixed == 0) {
            first_fix = cur_generation;
            printf("ff gen %d c %d\n",cur_generation,c);
        }
        if(cur_generation > generations / 2 && first_half_fixed < 0) {
            first_half_fixed = c_fixed;
        }
        c_fixed = c;
    }
}

void show_status(long cur_generation)
{
    printf("Gen %d: pop %d min %d fixed %d\n",cur_generation,cur_pop,min_pop,c_fixed);
}

void show_final()
{
    printf("After %d generations: population = %d (min %d)\n",generations, cur_pop, min_pop);
    printf("Total fixed genes: %d\n",c_fixed);
    printf(" first fix: generation %d\n",first_fix);
    printf(" fixed in first half %d\n",first_half_fixed);
    if(c_fixed > 1) {
        printf(" rate (post first): %g per generation\n", (c_fixed-1)/(float)(generations-first_fix));
        printf(" rate (last half of run): %g\n",(c_fixed-first_half_fixed)/(float)(generations*0.5f));
    }
    // For fun
    long icur,igene;
    long long c_mut=0;
    for(icur=0;icur<cur_pop;icur++) {
        for(igene=0;igene<GENE_BYTES;igene++) {
            c_mut += bit_count[p_genes[icur*GENE_BYTES+igene]];
        }
    }
    printf("Total mutations across population: %lld, mutations per organism %lld\n",
        c_mut, c_mut / cur_pop);
}

int main(int argc, char **argv)
{
    // Parse arguments
    if(argc<2) {
        // Print usage and exit
        printf("Usage: %s m [n [g [l [s [d]]]]]\n",argv[0]);
        printf("  m = mutation rate (capped at 1.0) per birth\n");
        printf("  n = population size [default: 1000]\n");
        printf("  g = generations (needs to be at least 4*n to get close to steady state) [default: 40*n]\n");
        printf("  l = fraction of mutations that are lethal (0.01 = 1%) [default: 0.01]\n");
        printf("  s = fraction of mutations that are strongly negative (10% lower survival chance) [default: 0.02]\n");
        printf("  d = fraction of mutations that are mildly negative (1% lower survival chance) [default: 0.03]\n");
        return(-1);
    }
    printf("argc %d\n",argc);
    m_rate = atof(argv[1]);
    pop_size = 1000;
    neg_l = 0.01;
    neg_s = 0.02;
    neg_m = 0.03;
    if(m_rate > 1.0f) m_rate = 1.0f;
    if(argc > 2) pop_size = atol(argv[2]);
    if(argc > 3) generations = atol(argv[3]);
    else generations = 40*pop_size;
    if(argc > 4) neg_l = atof(argv[4]);
    if(argc > 5) neg_s = atof(argv[5]);
    if(argc > 6) neg_m = atof(argv[6]);

    // Check negatives
    if(neg_l < 0 || neg_s < 0 || neg_m < 0) {
        printf("Negative mutation probabilities cannot be less than 0. (%f %f %f)\n",
            neg_l,neg_s,neg_m);
        return(-2);
    }
    if(neg_l + neg_s + neg_m > 1.0f) {
        printf("Negative mutation probabilities cannot sum to more than 1.0: (%f %f %f)\n",
            neg_l,neg_s,neg_m);
        return(-3);
    }
    neg_lb = neg_l * GENE_BYTES;
    neg_sb = neg_s * GENE_BYTES;
    neg_mb = neg_m * GENE_BYTES;
    printf("neg_bytes %d %d %d\n",neg_lb,neg_sb,neg_mb);
    // Allocate memory for organisms
    p_genes = calloc(pop_size*2, GENE_BYTES);
    p_offspring = calloc(pop_size*2, GENE_BYTES);
    if(p_genes == NULL || p_offspring == NULL) {
        printf("Insufficient memory to host %d organisms and 2 x %d offspring, please try a lower population size.\n",
            pop_size);
        if(p_genes != NULL) free(p_genes);
        return(-4);
    }

    // Print our values before we start
    printf("Starting run with mutation rate = %f, population size = %d, running %d generations.\n",
        m_rate, pop_size, generations);
    printf(" Negative mutation rates: %f lethal, %f strongly negative, %f mildly negative\n",
        neg_l, neg_s, neg_m);

    // Seed RNG
    long tval = (long)time(NULL);
    printf(" srand(%d)\n",tval);
    srand(tval);


    // Set up the optimized bit counter array
    for(int i=0;i<=255;i++) {
        int bits = 0;
        int tmp = i;
        while(tmp != 0) {
            tmp &= (tmp-1);
            bits++;
        }
        bit_count[i] = bits;
    }

    // Initialize - p_genes already 0 from calloc
    cur_pop = pop_size;
    min_pop = pop_size;
    // Loop and evaluate
    for(long i=0;i<generations;i++) {
        do_breed();
        do_cull();
        if(i % 10 == 0) {
            do_evaluate(i);
        }
        if(i % 1000 == 0) {
            show_status(i);
        }
    }
    show_final();
    free(p_genes);
    free(p_offspring);

    // TODO: Wrap the above in a "# of experiments" loop later
}