mitosim.md

mitochondrial and y-chromosomal lineage simulator

Simulates a population over several generations, tracking the mitochondrial eve and y-chromosome adam of each individual.

Each generation every female has a chance to produce up to 12 children. Half are male and dropped from the total population. The stop condition for reproduction is experimentally determined to keep the total population stable. No bias is applied to any lineage.

The RNG is a linear shift register based on the highest period shift register listed on wikipedia.

#How to run:

created and tested in https://webassembly.studio/ . provide implementation for printNb to log the size of the population and how many unique lineages survive.

license

This work is licensed under a Creative Commons Attribution 4.0 International License.

simulation.c

#define LINE_NUM 15000
#define MAX_LINE_NUM (LINE_NUM*5)
int lines1[MAX_LINE_NUM];
int lines2[MAX_LINE_NUM];

int lines3[MAX_LINE_NUM];
int lines4[MAX_LINE_NUM];

unsigned int  seed = 0xcafebabe;

//LCG based on 
// https://en.wikipedia.org/wiki/Linear-feedback_shift_register#Example_polynomials_for_maximal_LFSRs
int rand(int bits){
  unsigned taps = 0xE10000;

  for(int i=0; i<bits;i++){
    unsigned int tmp = (seed & taps)>>16;
    seed = seed << 1;
    seed |= (tmp ^ (tmp >> 5) ^ (tmp >> 6) ^ (tmp >> 7))& 1; 
  }

  return seed & ((1<<bits)-1);
}

int randInt(int max){
  int res;
  int bitid = 0;

  if(max < 2)return 0;

  {
    int tmp = max;
    if(tmp>0x000ffff) {bitid+=16;tmp >>= 16;}
    if(tmp>0x00000ff) {bitid+=8; tmp >>=  8;}
    if(tmp>0x000000f) {bitid+=4; tmp >>=  4;}
    if(tmp>0x0000003) {bitid+=2; tmp >>=  2;}
    if(tmp>0x0000001) {bitid+=1; tmp >>=  1;}
  }

  do{
    res = 0;
    int mask = bitid;
    while(mask>=0){
      res|= rand(1)<<mask;
      mask--;
      if(res >= max){
        res = 0;
        mask = bitid;
      }
      
    }
  }while(res>=max);
  return res;
}

int countUniqueUnsorted(int* lines, int alive, int* scratch){
  for(int i = 0; i<LINE_NUM; i++){
    scratch[i] = 0;
  }
  for(int i = 0; i< alive; i++){
    scratch[lines[i]]++;
  }
  int unique = 0;
  for(int i = 0; i<LINE_NUM; i++){
    if(scratch[i] != 0)unique++;
  }return unique;
}


void shuffle(int* arr, int size){
  for(int i = 0; i< size-1; i++){
    int randidx = i + randInt(size-i);
    int tmp = arr[i];
    arr[i] = arr[randidx];
    arr[randidx] = tmp;
  }
}

extern void printNb(int, int, int, int);

WASM_EXPORT
void simulate() {
  for(int i = 0; i<LINE_NUM; i++){
    lines1[i] = i+1;
    lines3[i] = i+1;
  }
    
  
  int alive_girls = LINE_NUM;
  int* alive_girl_lines = lines1;
  int* dest_girl_lines = lines2;

  int alive_boys = LINE_NUM;
  int* alive_boy_lines = lines3;
  int* dest_boy_lines = lines4;

  
  int generations = 1;
  for(; generations<1000; generations++){
    int dest_alive_girls = 0;
    int dest_alive_boys = 0;
    //randomly pair up everyone
    shuffle(alive_girl_lines, alive_girls);
    shuffle(alive_boy_lines, alive_boys);
    int alive = alive_girls;
    if(alive_boys<alive)alive = alive_boys;

    for(int i = 0; i<alive; i++){
      for(int j = 0; j<12; j++){
        int bits = 6;
        int cutoff = 21;
        if(alive>LINE_NUM){
          //population too high; lower amount of children per pair (slight hack to keep population stable)
          if(((rand(bits)) < cutoff+1))break;
        }else {
          if(((rand(bits)) < cutoff))break;
        }
        if(rand(1)==0)
          dest_girl_lines[dest_alive_girls++]=alive_girl_lines[i];
        else 
          dest_boy_lines[dest_alive_boys++] = alive_boy_lines[i];
        if(dest_alive_girls == MAX_LINE_NUM)break;//can't hold any more children
        if(dest_alive_boys == MAX_LINE_NUM)break;//can't hold any more children
      }
      if(dest_alive_girls == MAX_LINE_NUM)break;//no point
      if(dest_alive_boys == MAX_LINE_NUM)break;//no point
    }
    alive_girls = dest_alive_girls;
    alive_boys = dest_alive_boys;

    int* tmp = alive_girl_lines;
    alive_girl_lines = dest_girl_lines;
    dest_girl_lines = tmp;

    tmp = alive_boy_lines;
    alive_boy_lines = dest_boy_lines;
    dest_boy_lines = tmp;

    
    int unique_girls = countUniqueUnsorted(alive_girl_lines, alive_girls, dest_girl_lines);
    int unique_boys = countUniqueUnsorted(alive_boy_lines, alive_boys, dest_boy_lines);
    
    printNb(unique_girls, alive_girls, unique_boys, alive_boys);
    if(unique_girls == 1 && unique_boys == 1)break;
  }

}