shiba.c

Go to the documentation of this file.

#include "shiba.h"

void shiba(phylo phy)
{
  // pthreads preparation:
  pthread_t *thread = (pthread_t *) malloc(RUN_BATCH * sizeof(pthread_t));
  pthread_attr_t attr;
  void *status; // for pthread_join
  pthread_attr_init(&attr);
  pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);

  // these are global variables passed between functions. 
  // Shrink the scope if poss, and rename
  record = mem2d_i(Lineages, Spaces);
  success = 0;
  topresent = 0;
  maxarea = findMaxArea();
  maxdist = findMaxDist();

  // set disp, surv functions
  PSurv = mem2d_d(Times, Spaces);
  PDisp = mem3d_d(Times, Spaces, Spaces);
  for (int t = 0; t < Times; t++)
    for (int s = 0; s < Spaces; s++) {
      PSurv[t][s] = pSurvival(Area[t][s] /  maxarea);
      for (int s2 = 0; s2 < Spaces; s2++) 
        PDisp[t][s][s2] = pDispersal(Dist[t][s][s2] /  maxdist);
    }
  
  // initialize srandom with time
  srandom(time(NULL));

  printf("\n# SHIBA (Simulated Historical Island Biogeography Analysis)\n\n");
  
  long int batch;
  // run in batches of RUN_BATCH 
  for (batch = 0; batch < (long int)(Cfg.maxRuns/RUN_BATCH);batch++) {
    // start threads
    for (long int i = 0; i < RUN_BATCH ; i++)
      pthread_create(&thread[i], &attr, biogeo, NULL);

    // collect:
    for (long int i = 0; i < RUN_BATCH ; i++)
      pthread_join(thread[i], &status);

    // Reporting (and adjustment)
    fprintf(stderr,"\r> runs: %9ld;    to present: %8ld;    successes: %5ld", 
            batch * RUN_BATCH , topresent , success);
    fflush(NULL);  // Flush the buffer

    // check for successes
    if (success >= Cfg.stopAfterSuccess) break;
  }

  // Final reporting:
  fprintf(stderr, "\n\n"); fflush(NULL);
  printf("## Successes : %ld (in %ld runs)\n\n", success, batch * RUN_BATCH );
  printSuccessAll(phy);

  // free3d_i(locn , Lineages, Times);
  free2d_i(record, Lineages);

  free(thread); free(status);
  pthread_attr_destroy(&attr);
  pthread_mutex_destroy(&mymutex);

}

double findMaxArea()
{
  // finds total area on the landscape
  // was - max island size.  Now unlikely that any area get 100% survival
  double tot, max = 0.0;
  for (int t = 0; t < Times; t++) {
    tot = 0;
    for (int a = 0; a < Spaces; a++) tot += Area[t][a];
    if (tot > max) max = tot;
  }
  return max;
}

double findMaxDist()
{
  double max = 0.0;
  for (int t = 0; t < Times; t++) 
    for (int a = 0; a < Spaces; a++) 
      for (int b = 0; b < Spaces; b++) 
        if (Dist[t][a][b] > max) max = Dist[t][a][b];
  return max;
}

void *biogeo()
{

  int ***locn = mem3d_i(Lineages, Times, Spaces);



  int ss = 0;
  while (ss < Cfg.nStartSpaces) {
    int rnd = (int)((double) Spaces * (double) random() /
                    (double)(RAND_MAX+1.0));


    if ((!locn[0][0][rnd]) && 
        (Area[0][rnd] > 0.0) && 
        (Cfg.startSpace[rnd])) {
      locn[0][0][rnd] = 1;
      ss++;
    }
  }

  for (int t = 0; t < Times-1; t++)  { 
    // print array at start
    if (Cfg.verbose==1) 
      for (int i = 0; i < Lineages; i++)
        if (LineagePeriod[i][t]) {
          printf("t%d l%d : ", t, i);
          for (int x = 0; x < Spaces; x++) printf("%d", locn[i][t][x]);
          printf("\n");
        }
  
    for (int l = 0; l < Lineages; l++) {

      if (LineagePeriod[l][t]) {

        int linck = 0; // the number of surv lineages from one period to next
        
        for (int s = 0; s < Spaces; s++) {


          if (locn[l][t][s]) {

            for (int i = 0; i < Spaces ; i++) {
              if ( ( ((double) random() / (double)(RAND_MAX+1.0))
                     < PDisp[t][s][i] ) && (i != s) && (Area[t][i] > 0.0)) {
                locn[l][t+1][i] = 1;
                
                if (Cfg.verbose==1) printf("               l%d disperses from s%d to s%d\n",l,s,i);
                linck++;
              }
            }
            
            if (((double) random() / (double)(RAND_MAX+1.0)) < PSurv[t][s]) {
              locn[l][t+1][s] = 1;
              linck++;
            }
            else if (Cfg.verbose==1) printf(
              "               l%d goes locally extinct in s%d\n", l , s);
            
          }
        }

        if (!linck) {
          if (Cfg.verbose==1) printf("               l%d has died out\n",l);
          if (Cfg.verbose==1) printf("----------\n");
          free3d_i(locn, Lineages, Times);
          pthread_exit(NULL);
          //return;
        }


        if (!LineagePeriod[l][t+1]) {


          int places = 0;      //number of spaces occupied by daughters of `l`
          int placemapping[Spaces]; // mapping from 0..places to 'real' spaces
          for (int s = 0; s < Spaces; s++)
            if (locn[l][t+1][s]) {
              placemapping[places] = s;  
              places++;
            }
          int placeck[places]; // has the place got a lineage in it? 
          for (int h = 0; h < places; h++) placeck[h] = 0;
          int fullplaces = 0;       // number of places filled
          
          
          for (int q = 0; q < LineageDaughtersN[l]; q++) {
            int placed = 0; // has this daughter been placed?
            while (!placed) {
              
              int rnd = (int) ( (double) places * (double) random() / 
                                (double) (RAND_MAX+1.0));
              if (placeck[rnd]) {
                if (fullplaces >= places) {
                  locn[LineageDaughters[l][q]][t+1][placemapping[rnd]] = 1;
                  placed = 1;
                  if (Cfg.verbose==1) printf(
                    "               l%d ends; l%d begins in s%d\n", l, 
                    LineageDaughters[l][q], placemapping[rnd]);
                }
              } else {
                locn[LineageDaughters[l][q]][t+1][placemapping[rnd]] = 1;
                placeck[rnd] = 1;
                placed = 1;
                fullplaces++;
                if (Cfg.verbose==1) printf(
                  "               l%d ends; l%d begins in s%d\n", 
                  l, LineageDaughters[l][q], placemapping[rnd]);

              }
            }
          }
        }
      }
    }
  }

  // final layout
  if (Cfg.verbose==1)
    {
      for (int i = 0; i < Lineages; i++)
        if (LineagePeriod[i][Times-1]) {
          printf("t%d l%d : ", Times-1, i);
          for (int x = 0; x < Spaces; x++) printf("%d", locn[i][Times-1][x]);
          printf("\n");
        }

      //printExtant(Lineages);
    }

  
  int diffs = 0;  // the count of differences between sim and obs
  for (int l = 0; l < Lineages; l++) {
    if (LineagePeriod[l][Times-1]) {
      for (int s = 0; s < Spaces; s++) {
        if (locn[l][Times-1][s] != LineageExtant[l][Times-1][s]) diffs++;
      }
    }
  }
  

  // lock mutex
  pthread_mutex_lock(&mymutex);
  topresent++;

  if (!diffs) {
    success++;
    for (int l = 0; l < Lineages; l++) {
      for (int t = 0; t < Times; t++) {
        if (LineagePeriod[l][t]) {
          for (int s = 0; s < Spaces; s++)
            record[l][s] += locn[l][t][s];
          break;
        }
      }
    }
  }
  pthread_mutex_unlock(&mymutex);

  if (Cfg.verbose==1) printf("-----------\n");

  free3d_i(locn, Lineages, Times);
  pthread_exit(NULL);
}


void printSuccessAll(phylo p)
{
  printf("## Number of runs in which the base of each lineage was in each space:\n   (*= extant distrib)\n\n");
  // printf("lin  bor   wsu   ngu   phi   luz\n");
  printf("    Space       |  ");
  for (int k = 0; k < Spaces; k++) printf(" %3d ", k);
  printf("\n");
  printf("    ------------+--");
  for (int k = 0; k < Spaces; k++) printf("-----");
  printf("\n");
      
  for (int i = 0; i < Lineages; i++) {
    printf("    Lineage %3d |  ", i);
    for (int k = 0; k < Spaces; k++) {
      if (LineageExtant[i][Times-1][k])
        printf(" %3d*", record[i][k]);
      else  printf(" %3d ", record[i][k]);
    }
    printf("  --> ");
    if(p.taxon[i]) printf("%s", p.taxon[i]);
    for(int l = 0; l < LineageDaughtersN[i]; l++)
      printf("l%d ", LineageDaughters[i][l]);
    printf("\n");
  }
  printf("\n");
}

double pDispersal(double x) {
  return Cfg.probDispA * powf( 10.0, -1.0 * Cfg.probDispB * x ) ;
}

double pSurvival(double x) {  
  return Cfg.probSurvA * (log10f(( x * ( powf(10, Cfg.probSurvB) -1.0 )) +1.0 ) 
                           / Cfg.probSurvB );
}