/* contevol.c, written by Joe Felsenstein
 * (c) Copyright 1999 by Joe Felsenstein.
 * Simulates a quantitative character under natural selection
 * version of 21 December 1999
 * Permission is granted to copy and use this program provided this
 * notice (including copyright and program authors) is not removed
 * and no fee is charged for the program.
 */

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


/*function prototypes for ANSI C compatibility */

#ifndef OLDC
void uppercase(char *);
double dgetnum();
long lgetnum();  
long lgetnum();
double randreal();
char get_parameters(int, long *, long *, double *,
                     double *, double *, double *, long *);
void make_fitnesses(long, double, double, double);
void do_n_generations(int, long);
int main(int, char **);
#endif

#define MAX_POP_SIZE    10000
#define LOCI            5
#define LOCI2           11
#define MAX_GEN_RUN     10000
#define SCREENLINES     24

#define TRUE  1
#define FALSE 0
typedef char Boolean;

#define LARGE_BUF_LENGTH 500
char Buffer[LARGE_BUF_LENGTH];

/*
 * Global variables  :(
*/

typedef int DisplayType;

DisplayType disp_type;
typedef long genotype_array[LOCI];

genotype_array x[MAX_POP_SIZE]; /* the array individuals in a generation */
genotype_array y[MAX_POP_SIZE]; /* the array individuals in next generation */
long phen[MAX_POP_SIZE];        /* the phenotypes in the next generation */
long xx[LOCI];                  /* temporary array for an individual */
double w[LOCI2];                /* fitnesses of phenotypes */
long hist[LOCI2];               /* histogram of phenotypes */
long t1;                        /* the beginning generation for a run */
long npop;                      /* the size of the population */
long init_pop_size;             /* initial population size */
long g9;                        /* generations to run (just 1 for now )*/
long peaks;                     /* number of peaks of fitness curve */
Boolean startover;              /* Boolean for restarting same case */
Boolean change;                 /* Boolean for changing the parameters */
Boolean menuagain;              /* Boolean for reconsidering fitnesses */


void uppercase(ch)
     char *ch;
{
  /* convert ch to upper case -- either ASCII or EBCDIC */
  *ch = isupper(*ch) ? *ch : toupper(*ch);
}  /* uppercase */


double dgetnum() 
{
  double num;

  fgets(Buffer, LARGE_BUF_LENGTH, stdin);
  while ( sscanf(Buffer, "%lg", &num) != 1) {
    printf("\nPlease enter a number: ");
    fgets(Buffer, LARGE_BUF_LENGTH, stdin);
  }
  
  return num;
} /* dgetnum */


long lgetnum() 
{
  long num;

  fgets(Buffer, LARGE_BUF_LENGTH, stdin);
  while ( sscanf(Buffer, "%ld", &num) != 1) {
    printf("\nPlease enter a number: ");
    fgets(Buffer, LARGE_BUF_LENGTH, stdin);
  }

  return num;
} /* lgetnum */


double randreal()
{ /* the random number generator */
#ifdef __MWERKS__
  return (((double)(rand())) / (double) RAND_MAX);
#else
#ifdef __TURBOC__
  return (((double)(rand())) / (double) RAND_MAX);
#else
#ifdef __WATCOMC__
  return (((double)(rand())) / (double) RAND_MAX);
#else  
  return (((double)(random())) / (double)2147483647); /* div by RAND_MAX */
#endif
#endif
#endif
} /* randreal */


void cleer()
{ /* Clear screen */
#ifndef __MWERKS__
   printf("\033[2J\033[H");
#endif
} /* cleer */


Boolean get_parameters(disp_type, pop_size, peaks, opt1, opt2, sel, init_phen,
                       gen_to_run)
     DisplayType disp_type;
     long *pop_size;
     long *peaks;
     double *opt1;
     double *opt2;
     double *sel;
     double *init_phen;
     long *gen_to_run;

{ /* read in values of parameters for run */
  char answer;
  static long prev_gen_to_run = -1;
  Boolean done;

  if (prev_gen_to_run == -1) {
    prev_gen_to_run = *gen_to_run;
  }
  else {
    *gen_to_run = prev_gen_to_run;
  }

  done = FALSE;
  while (!done) { 
    cleer();
    printf("\nContevol settings\n");
    printf("Option Character       What it is         Its current value\n");
    printf("   P     Number of peaks in fitness curve  %ld\n", *peaks);
    if (*peaks == 1)
      printf("   O              Optimum character value  %lf\n", *opt1);
    else {
      printf("   1           Optimum character value #1  %lf\n", *opt1);
      printf("   2           Optimum character value #2  %lf\n", *opt2);
    }
    printf("   S                Strength of selection  %lf\n", *sel);
    printf("   N                     Population size:  %ld\n", *pop_size);
    printf("   I           Initial mean of character:  %lf\n", *init_phen);
    printf("   G         Generations to next display:  %ld\n", *gen_to_run);
   printf("Are these settings correct? ([Y], Q (quit), or option to change): ");
    
    fgets(Buffer, LARGE_BUF_LENGTH, stdin);
    answer = Buffer[0];
    
    uppercase(&answer);
    switch (answer) {
    case 'Y': case '\n': case '\0': /* done with setting the parameters */
      t1 = 0;                        /* reset the beggining time */
      return FALSE; /* we are NOT done with the run */
      break;
    case 'Q': /* done with the simulation */
      done = TRUE;
      return TRUE;  /* the simulation is over */
      break;
    case 'P': /* Number of peaks */
      printf("Number of peaks in fitness curve (1 or 2)? ", *peaks);
      *peaks = lgetnum();
      break;
    case 'O': /* Optimum phenotype */
      printf("Optimum phenotype? ", *opt1);
      *opt1 = dgetnum();
      break;
    case '1': /* Optimum phenotype #1 */
      if (*peaks == 1)
        printf("Optimum phenotype? ", *opt1);
      else
        printf("Optimum phenotype #1? ", *opt1);
      *opt1 = dgetnum();
      break;
    case '2': /* Optimum phenotype #1 */
      if (*peaks == 1)
        printf("Optimum phenotype? ", *opt1);
      else
        printf("Optimum phenotype #2? ", *opt2);
      *opt2 = dgetnum();
      break;
    case 'S': /* Selection strength */
      printf("Selection strength (big is weak)? ", *sel);
      *sel = dgetnum();
      break;
    case 'N': /* population */
      printf("\nNew population size (was %d)? ", *pop_size);
      *pop_size = lgetnum();
      while (((unsigned long)*pop_size > MAX_POP_SIZE) && 
             ((unsigned long)*pop_size <= 0)) {
        printf("\nPlease enter a number between %d and %d: ", 1, MAX_POP_SIZE);
        *pop_size = lgetnum();
      }
      init_pop_size = *pop_size;
      break;
    case 'I': /* initial phenotype mean */
      printf("Initial mean character value? ", *init_phen);
      *init_phen = dgetnum();
      while ((*init_phen < 0.0) && (*init_phen > (double)LOCI2)) {
        printf("\nPlease enter a number between 0.0 and %lf: ",
                                                   (double)(2*LOCI));
        *init_phen = dgetnum();
      }
      break;
    case 'G': /* generations to run */
      printf("How many generations until next histogram is displayed: ");
      *gen_to_run = lgetnum();
      while (((unsigned long)*gen_to_run > MAX_GEN_RUN) ||
             ((unsigned long)*gen_to_run < 0)) {
        printf("\nPlease enter a number between 0 and %d: ", MAX_GEN_RUN);
        *gen_to_run = lgetnum();
      }
      prev_gen_to_run = *gen_to_run;
      break;
    default:
      printf("\nNot an option\n");
      break;
    }
  } /* end loop */
} /* get_parameters */


void make_fitnesses(peaks, opt1, opt2, sel)
   long peaks;
   double opt1, opt2, sel;
{ /* make a table of fitnesses of phenotypes */
  long i, j, k, l, m;
  double maxw, ww, f;
  char c;

  for (i=0; i <= 2*LOCI; i++) {     /* will be a one or two-peaked Gaussian */
    w[i] = exp(-0.5*(((double)i-opt1)*((double)i-opt1)/sel));
    if (peaks == 2)
      w[i] += exp(-0.5*(((double)i-opt2)*((double)i-opt2)/sel));
  }
  maxw = 0.0;                         /* find maximum fitness */
  for (i=0; i <= 2*LOCI; i++)
    if (w[i] > maxw) maxw = w[i];
  for (i=0; i <= 2*LOCI; i++)         /* divide all fitnesses by it */
    w[i] = w[i] / maxw;
  cleer();
  for (i = SCREENLINES-2; i >= 1; i--) {        /* go down screen */
    if ((SCREENLINES-2-i) % (long)((SCREENLINES-2)/10) == 0)
                                      /* print fitness scale at left */
      printf("%4.2lf", (double)(i)/(double)(SCREENLINES-2));
    else
      printf("    "); 
    printf(" |   ");                  /* print left edge */
    m = (LOCI+1)/2;
    for (j = m; j <= 12*LOCI+m; j++) {
      k = (j-m) / (LOCI+1);          /* which fitness value is to left */
      l = (j-m) % (LOCI+1);          /* how many spaces left of here */
      f = (double)l / (LOCI+1);      /* weight on fitness to right, 0.0-1 */
      if (l == 0)                    /* if we're at a phenotype value ... */
        ww = w[k];                   /* get its fitness */
      else                           /* if we're in between ... */
        ww = (1-f)*w[k] + f*w[k+1];  /* interpolate, */
      if ((long)(ww*(SCREENLINES-2)+0.5) == i)  {
        if (l == 0)
          printf("O");
        else
          printf(".");
      }
      else
        printf(" ");
    }
    printf("\n");
  }
  printf("     L-");                             /* print bottom row */
  for (i = 0; i <= 2*LOCI; i++)             /* print bottom row */
    if (i < 10)
       printf("--%1ld---", i);
    else
       printf("-%2ld---", i);
  printf("\n");                             /* print bottom row */
  printf(
   "to continue, press Enter key, to stop press Q, to do menu again press M ");
  fgets(Buffer, LARGE_BUF_LENGTH, stdin);  /* wait till user presses Enter */
  c = Buffer[0];
  menuagain = (toupper(c) == 'M');
  if (toupper(c) == 'Q') {
    printf("\nThat's all, folks!\n\n");
    exit(0);
  }
} /* make_fitnesses */


void do_n_generations(disp_type, g9)
     DisplayType disp_type;
     long g9;                       /* generations to run (actually 1) */

{
  long i, j, t, k, l, pphen, maxhist;
  double ww;
  long where[LOCI2];
  char c;

   /* The section below is where all the quantitative genetics happens.
      If you are trying to understand the simulation this loop is all
      you need to know about.  The rest is input/output */

  for (t = 0; t < g9; t++) {
    for (i = 0; i <= 2*LOCI; i++)         /* intialize histogram */
      hist[i] = 0;
    for (i = 1; i <= npop; i++) {         /* for each resulting individual */
      do {
        pphen = 0;
        for (j = 1; j <= LOCI; j++) {     /* segregate a locus from 2 parents */
          xx[j] = (long)(randreal()+(double)x[(long)(randreal()*npop)+1][j]/2)
                + (long)(randreal()+(double)x[(long)(randreal()*npop)+1][j]/2);
          pphen = pphen + xx[j];          /* add up the phenotype */
        }
        ww = w[pphen];
      } while (randreal() > ww);          /* until an offsprint survives */
      for (j=1; j <= LOCI; j++)            /* put genotype in array y */
        y[i][j] = xx[j];
    }                                     /* end of loop over individuals */
    for (i = 1; i <= npop; i++) {         /* copy y back into x */
      phen[i] = 0;
      for (j = 1; j <= LOCI; j++) {
        x[i][j] = y[i][j];
        phen[i] += y[i][j];               /* and (re-)make the phenotypes */
      }
      hist[phen[i]]++;                    /* and add to histogram */
    }

     /* end of the important quantitative genetics stuff */

  }  /* end of for t < g9 -- end of running through the (1) generations */
  maxhist = 0;                             /* find biggest # in histogram */
  for (i = 0; i <= 2*LOCI; i++)
    if (hist[i] > maxhist) maxhist = hist[i];
  cleer();
  for (i = SCREENLINES/2-2; i >= 1; i--) {        /* go down screen */
    for (k = 1; k <= 2; k++) {             /* two rows per layer of histogram */
      printf(" |");                        /* print left edge */
      for (j = 0; j <= 2*LOCI; j++)        /* intialize locations of indiv's */
         where[j] = 0;
      for (j = 0; j <= 2*LOCI; j++) {
        if (k == 1)
          printf(" ");
        else
          printf("_");
        if ((long)(((double)hist[j]/maxhist)*(SCREENLINES/2-2)+0.5) >= i)  {
                                           /* if there's someone here */
          do {
            where[j]++;                    /* check next one */
          } while (phen[where[j]] != j);   /* loop until find one */
          for (l = 1; l <= LOCI; l++) {    /* prints genes */
            if ((x[where[j]][l] == 2) || ((x[where[j]][l] == 1) && (k == 1)))
              printf("%c", (char)(64+l));
            else
              printf("%c", (char)(96+l));
          }
        }
        else {
          for (l = 1; l <= LOCI; l++)       /* no individual here */
            printf(" ");
        }
      }
      printf("\n");
    }
  }
  printf(" L-");                             /* print bottom row */
  for (i = 0; i <= 2*LOCI; i++)        
    if (i < 10)
       printf("--%1ld---", i);
    else
       printf("-%2ld---", i);
  printf("\n");
  printf("  ");                             /* print bottom row */
  for (i = 0; i <= 2*LOCI; i++)             /* print bottom row */
    printf(" %3ld  ", hist[i]);
  printf("\n");
  printf("to continue, press Enter key, to stop press Q  ");
  fgets(Buffer, LARGE_BUF_LENGTH, stdin);  /* wait till user presses Enter */
  c = Buffer[0];
  if (toupper(c) == 'Q') {
    printf("\nThat's all, folks!\n\n");
    exit(0);
  }
  if (toupper(c) == 'S')
    startover = TRUE;
  else 
    startover = FALSE;
  if (toupper(c) == 'C')
    change = TRUE;
  else 
    change = FALSE;
} /* do_n_generations */


int main(argc, argv)
     int argc;
     char *argv[];
{  /* main */
  long i, j;
  long g9;
  double p0, opt1, opt2, sel, init_phen;

  Boolean cont = TRUE;
  Boolean done = FALSE;

  /* default values for the parameters */
  peaks = 1;   /* no. of peaks of fitness curve */
  npop  = 100;  /* population */
  p0 = 0.5;  /* initial freq. of A */
  g9 = 1;    /* generations to run */
  opt1 = 5.0;  /* optimum value #1 */
  opt2 = 5.0;  /* optimum value #2 */
  sel = 10.0;  /* (inverse) strength of selection */
  init_phen = 5.0;      /* initial mean phenotype */
  startover = FALSE;    /* boolean for restarting same case */
  change = FALSE;       /* Boolean for getting new parameters */
  menuagain = FALSE;    /* Boolean for re-doing fitnesses */

  printf("Contevol  (c) Copyright Joe Felsenstein 1999\n");

#ifdef __MWERKS__
  srand((int)time(NULL));
#endif
#ifndef __MWERKS__
  srandom((int)time(NULL));
#endif

  disp_type = 1;

  do {
    if (!startover) {
      do {
        done = get_parameters(disp_type, &npop, &peaks, &opt1, &opt2,
                                       &sel, &init_phen, &g9);
        if (!done)
          make_fitnesses(peaks, opt1, opt2, sel);
      } while (menuagain);
    }
    if (!done) {           /* initialize the starting population */
      p0 = init_phen/(2*LOCI);
      for (i = 1; i <= npop; i++)
        for (j = 1; j <= LOCI; j++) {
           x[i][j] = (long)(p0+randreal()) + (long)(p0+randreal());
        }
      cont = TRUE;
      while (cont) {
        do_n_generations(disp_type, g9);
        /* update for next run and prompt to run more or stop */
        t1 += g9;
        if (g9 && !startover && !change)
          cont = TRUE;
        else {
          cont = FALSE;
        }
      }
    }
  } while (!done);
  printf("\nThat's all, folks!\n\n");
  exit(0);
}