/* sex2lo.c, May 5, 1989, Comparison of Asexual and Sexual Populations, jr */
/* Two-state Markov environment for two loci with two alleles per locus.   */

/* Selection strength is read from stdin. Provided the strength is > 0,    */
/* recombination, transition probs., mutation, run len., and num. of       */
/* trials are readin, then the ensemble of sexual/asexual comparisons is   */
/* computed. For each ensemble the expectation and standard error of       */
/* the ratio of the sexual to asexual population's arithmetic mean wbar,   */
/* standard deviation of wbar, and geometric mean of wbar are determined.  */
/* The population configurations at the end of each run together with      */
/* the ensemble statistics are printed.  Ensembles are reinitialized       */
/* and produced for each value of selection read from stdin until a        */
/* selection strength of 0 is encountered.  So, make a file of desired     */
/* parameters, six per line: selection, r, p12, p21, mutation, runlen,     */
/* and trials. Call this file "sex.in", and then usage is:                 */
/*                  $ sex2lo <sex2lo.in >sex2lo.out                        */
/* Ensemble summaries extracted with:                                      */
/*           $ grep "*" <sex2lo.out | tr "*" " " >sex2lo.sum               */
/* Sample parameter input file, sex2lo.in:                                 */
/*                    2 0.5 0.5 0.5 1.0e-8 1000 50                         */
/*                    0                                                    */
/* If r = 0.5, then the two loci are on different chromosomes.             */
/* If r = 0.0, then the two loci are adjacent on the same chromosome.      */
/* If p12 + p21 = 1 then the serial correlation between states is zero.    */
/* If p12 = p21 then both states are equally likely.                       */
/* In general: p2 = p12/(p12 + p21), stationary probability of bad env.    */
/*             p1 = p21/(p12 + p21), stationary prob. of good environment  */
/*            rho = 1 - p12 - p21, serial correlation between env. states  */

/* To compile without libblas.a: cc -O -f2 -o sex2lo sex2lo.c -lfm -lfc    */
/* With libblas.a: cc -O -f2 -DBLAS -o sex2lo sex2lo.c -lblas -lfm -lfc    */

/* genotype ordering ----                                                  */
/*           0: AB/AB, 1: AB/Ab, 2:AB/aB, 3:AB/ab, 4:Ab/Ab,                */
/*           5: Ab/aB, 6: Ab/ab, 7:aB/aB, 8:aB/ab, 9:ab/ab.                */

#define X11INITIAL 0.0        /* A2 and B2 alleles initially fixed in both */
#define X12INITIAL 0.0                   /* sexual and asexual populations */
#define X13INITIAL 0.0
#define X14INITIAL 0.0
#define X22INITIAL 0.0
#define X23INITIAL 0.0
#define X24INITIAL 0.0
#define X33INITIAL 0.0
#define X34INITIAL 0.0
#define X44INITIAL 1.0           /* Initial abundance for both populations */
#define SIZEINITIAL 0.5   /* equals half of what the universe will support */
#define FITNESS11 1   /* initial fitness state is "good" for all genotypes */
#define FITNESS12 1
#define FITNESS13 1
#define FITNESS14 1
#define FITNESS22 1
#define FITNESS23 1
#define FITNESS24 1
#define FITNESS33 1
#define FITNESS34 1
#define FITNESS44 1

                                                       /* system resources */
#include <stdio.h>           /* header file of standard system definitions */
extern double drand48();  /* returns double uniformly distributed on [0,1) */
extern void srand48();    /* seeds random number generator, call with SEED */
#define SEED 0L        /* seed random number generator with a long integer */
extern double sqrt();   /* sqrt(x) returns square root of x; x is a double */
extern double pow();   /* pow(x,y) returns x^y, where x,y are double > 0.0 */
extern double exp();            /* exp(x) returns e^x; where x is a double */
extern double log();     /* log(x) returns natural log of x; x is a double */
extern void printf();  /* standard function for formatted output to stdout */
extern void scanf();   /* standard function for formatted input from stdin */

                                              /* linear algebra procedures */
#ifdef BLAS                  /* define if libblas.a available on IBM RT PC */
extern double ddot_();                                      /* dot product */
extern void dscal_();                             /* scalar multiplication */
extern void dcopy_();                                       /* vector copy */
#else                         /* otherwise simulate the libblas.a routines */

double ddot_(p_n, p_x, p_incx, p_y, p_incy)                 /* dot product */
     int *p_n, *p_incx, *p_incy;     /* returns double precision x[] . y[] */
     double *p_x, *p_y;
     {
     register int i, n;
     register double *x, *y;
     register double ddot = 0;
     x = p_x;
     y = p_y;
     n = *p_n;
     for (i = 0; i < n; ++i)
          {
          ddot += (*x++)*(*y++);
          }
     return ddot;
     }

void dscal_(p_n, p_scaler, p_x, p_incx)
     int *p_n, *p_incx;                           /* scalar multiplication */
     double *p_scaler, *p_x;               /* replaces x[] with scaler*x[] */
     {
     register int i, n;
     register double *scaler, *x;
     scaler = p_scaler;
     x = p_x;
     n = *p_n;
     for (i = 0; i < n; ++i)
          {
          *x = (*scaler)*(*x);
          ++x;
          }
     }

void dcopy_(p_n, p_x, p_incx, p_y, p_incy)                  /* vector copy */
     int *p_n, *p_incx, *p_incy;                         /* sets y[] = x[] */
     double *p_x, *p_y;
     {
     register i, n;
     register double *x, *y;
     x = p_x;
     y = p_y;
     n = *p_n;
     for (i = 0; i < n; ++i)
          {
          (*y++) = (*x++);
          }
     }
#endif

int dimen = 10;                 /* dimension needed for libblas.a routines */
int inc = 1;              /* index increment needed for libblas.a routines */

void dschur_(p_n, p_x, p_incx, p_y, p_incy)               /* Schur product */
     int *p_n, *p_incx, *p_incy;            /* replaces y[] with x[] o y[] */
     double *p_x, *p_y;
     {
     register int i, n;
     register double *x, *y;
     x = p_x;
     y = p_y;
     n = *p_n;
     for (i = 0; i < n; ++i)
          {
          *y = (*x++)*(*y);
          ++y;
          }
     }

void dmxcv_(p_n, p_a, p_inca, p_x, p_incx, p_y, p_incy)   /* square matrix */
     int *p_n, *p_inca, *p_incx, *p_incy;           /* times column vector */
     double *p_a, *p_x, *p_y;             /* replaces y[] with a[][] * x[] */
     {
     register int i, n;
     register double *row_a, *x, *y;
     row_a = p_a;
     x = p_x;
     y = p_y;
     n = *p_n;
     for (i = 0; i < n; ++i)
          {
          *y++ = ddot_(p_n,row_a,p_inca,x,p_incx);
          row_a += n;
          }
     }

               /* variables used for each run, reinitialized by init_run() */
int fitrecord[10] = {0,};                /* record of previous environment */
double fitness[10] = {0.0,};  /* genotype fitnesses in current environment */
double sex_genotypes[10] = {0.0,};    /* genotype frequencies, sexual pop. */
double asex_genotypes[10] = {0.0,};  /* genotype frequencies, asexual pop. */
double sex_size = 0.0;                           /* sexual population size */
double asex_size = 0.0;                         /* asexual population size */
double sex_wbar = 0.0;                 /* sexual population's mean fitness */
double asex_wbar = 0.0;               /* asexual population's mean fitness */
double sex_amwbar = 0.0;               /* arithmetic mean of sexual's wbar */
double asex_amwbar = 0.0;             /* arithmetic mean of asexual's wbar */
double sex_sdwbar = 0.0;            /* standard deviation of sexual's wbar */
double asex_sdwbar = 0.0;          /* standard deviation of asexual's wbar */
double sex_gmwbar = 0.0;                /* geometric mean of sexual's wbar */
double asex_gmwbar = 0.0;              /* geometric mean of asexual's wbar */
double ratio_amwbar = 0.0;     /* ratio at end of run of sexual to asexual */
double ratio_sdwbar = 0.0;       /* population's arithmetic mean, standard */
double ratio_gmwbar = 0.0;                /* deviation, and geometric mean */

   /* ensemble parameters and statistics, reinitialized in init_ensemble() */
double sercorr = 0.0;   /* serial correlation between environmental states */
double probbad = 0.0;        /* stationary probability of a bad generation */
double wgood = 0.0;                        /* fitness in a good generation */
double wbad = 0.0;                          /* fitness in a bad generation */
double mutation[10][10] = {0.0,};              /* recurrent mutation rates */
double examratio = 0.0;     /* expectation and standard error for ratio of */
double seamratio = 0.0;          /* sexual to asexual arithmetic mean wbar */
double exsdratio = 0.0;     /* expectation and standard error for ratio of */
double sesdratio = 0.0;    /* sexual to asexual standard deviation of wbar */
double exgmratio = 0.0;     /* expectation and standard error for ratio of */
double segmratio = 0.0;           /* sexual to asexual geometric mean wbar */

         /* ensemble parameters, reinitialized from stdin or by assignment */
int selection = 0;      /* strength of natural selection, input from stdin */
double r = 0.0;                              /* recombination between loci */
double p12 = 0.0;    /* transition probability from good generation to bad */
double p21 = 0.0;    /* transition probability from bad generation to good */
double mutrate = 0.0;              /* recurrent within-locus mutation rate */
int timesup = 0;                      /* length of each run in generations */
int trials = 0;       /* number of trials to comprise statistical ensemble */

/* matrix of transition probabilities, pij, for transition from i to j is: */
/*               (1 - p12) p12                                             */
/*               p21 (1 - p21)                                             */

                                         /* model for fitness fluctuations */
void get_fitness()                    /* replace fitness[] with new values */
     {
     int i;
     for (i = 0; i < 10; ++i)
          {
          if (fitrecord[i] == 1)               /* last generation was good */
               {
               if (drand48() < p12)              /* moves from good to bad */
                    {
                    fitrecord[i] = 2;
                    fitness[i] = wbad;
                    }
               else                                          /* stays good */
                    {
                    fitrecord[i] = 1;
                    fitness[i] = wgood;
                    }
               }
          else /* fitrecord[i] == 2 */          /* last generation was bad */
               {
               if (drand48() < p21)              /* moves from bad to good */
                    {
                    fitrecord[i] = 1;
                    fitness[i] = wgood;
                    }
               else                                           /* stays bad */
                    {
                    fitrecord[i] = 2;
                    fitness[i] = wbad;
                    }
               }
          }
     fitrecord[3] = fitrecord[5];   /* force fitness equality for coupling */
     fitness[3] = fitness[5];        /* and repulsion double heterozygotes */
     }                       /* involves one superfluous call to drand48() */

void do_asex()         /* model for one generation of asexual reproduction */
     {
     double tmp[10],inv_wbar;
     asex_wbar = ddot_(&dimen,asex_genotypes,&inc,fitness,&inc);
     inv_wbar = 1.0/asex_wbar;
     dschur_(&dimen,fitness,&inc,asex_genotypes,&inc);
     dscal_(&dimen,&inv_wbar,asex_genotypes,&inc);
     asex_size = asex_wbar*asex_size;
     dmxcv_(&dimen,mutation,&inc,asex_genotypes,&inc,tmp,&inc);
     dcopy_(&dimen,tmp,&inc,asex_genotypes,&inc);
     }

void do_sex()           /* model for one generation of sexual reproduction */
     {
     double tmp[10],inv_wbar,p1,p2,p3,p4;
     sex_wbar = ddot_(&dimen,sex_genotypes,&inc,fitness,&inc);
     inv_wbar = 1.0/sex_wbar;
     dschur_(&dimen,fitness,&inc,sex_genotypes,&inc);
     dscal_(&dimen,&inv_wbar,sex_genotypes,&inc);
     sex_size = sex_wbar*sex_size;
     dmxcv_(&dimen,mutation,&inc,sex_genotypes,&inc,tmp,&inc);
     dcopy_(&dimen,tmp,&inc,sex_genotypes,&inc);                /* meiosis */
     p1 = sex_genotypes[0] + 0.5*sex_genotypes[1] + 0.5*sex_genotypes[2]
          + 0.5*(1.0 - r)*sex_genotypes[3] + 0.5*r*sex_genotypes[5];
     p2 = 0.5*sex_genotypes[1] + sex_genotypes[4]
          + 0.5*(1.0 - r)*sex_genotypes[5] + 0.5*sex_genotypes[6]
          + 0.5*r*sex_genotypes[3];
     p3 = 0.5*sex_genotypes[2] + 0.5*(1.0 - r)*sex_genotypes[5]
          + sex_genotypes[7] + 0.5*sex_genotypes[8]
          + 0.5*r*sex_genotypes[3];
     p4 = 0.5*(1.0 - r)*sex_genotypes[3] + 0.5*sex_genotypes[6]
          + 0.5*sex_genotypes[8] + sex_genotypes[9] + 0.5*r*sex_genotypes[5];
     sex_genotypes[0] = p1*p1;                   /* random union of gametes */
     sex_genotypes[1] = 2*p1*p2;
     sex_genotypes[2] = 2*p1*p3;
     sex_genotypes[3] = 2*p1*p4;
     sex_genotypes[4] = p2*p2;
     sex_genotypes[5] = 2*p2*p3;
     sex_genotypes[6] = 2*p2*p4;
     sex_genotypes[7] = p3*p3;
     sex_genotypes[8] = 2*p3*p4;
     sex_genotypes[9] = p4*p4;
     }

void fill_universe()
     {  /* universe repopulated in proportion to last generation's success */
     double combined;
     combined = asex_size + sex_size;
     asex_size = asex_size/combined;
     sex_size = sex_size/combined;
     }

void init_run()
     {
     sex_genotypes[0] = X11INITIAL;        /* initial genotype frequencies */
     sex_genotypes[1] = X12INITIAL;
     sex_genotypes[2] = X13INITIAL;
     sex_genotypes[3] = X14INITIAL;
     sex_genotypes[4] = X22INITIAL;
     sex_genotypes[5] = X23INITIAL;
     sex_genotypes[6] = X24INITIAL;
     sex_genotypes[7] = X33INITIAL;
     sex_genotypes[8] = X34INITIAL;
     sex_genotypes[9] = X44INITIAL;
     sex_size = SIZEINITIAL;                    /* initial population size */
     asex_genotypes[0] = X11INITIAL;
     asex_genotypes[1] = X12INITIAL;
     asex_genotypes[2] = X13INITIAL;
     asex_genotypes[3] = X14INITIAL;
     asex_genotypes[4] = X22INITIAL;
     asex_genotypes[5] = X23INITIAL;
     asex_genotypes[6] = X24INITIAL;
     asex_genotypes[7] = X33INITIAL;
     asex_genotypes[8] = X34INITIAL;
     asex_genotypes[9] = X44INITIAL;
     asex_size = SIZEINITIAL;
     fitrecord[0] = FITNESS11;            /* environment when clock starts */
     fitrecord[1] = FITNESS12;
     fitrecord[2] = FITNESS13;
     fitrecord[3] = FITNESS14;
     fitrecord[4] = FITNESS22;
     fitrecord[5] = FITNESS23;
     fitrecord[6] = FITNESS24;
     fitrecord[7] = FITNESS33;
     fitrecord[8] = FITNESS34;
     fitrecord[9] = FITNESS44;
     sex_amwbar = asex_amwbar =                   /* within-run statistics */
     sex_sdwbar = asex_sdwbar = 
     sex_gmwbar = asex_gmwbar = 0.0;
     }

void run_sex_asex()    /* run with sexual and asexual populations together */
     {         
     int generation;

     for (generation = 0; generation < timesup; ++generation)
          {
          get_fitness();
          do_asex();
          do_sex();
          fill_universe();

          sex_amwbar += sex_wbar;                       /* sum of sex_wbar */
          asex_amwbar += asex_wbar;                    /* sum of asex_wbar */
          sex_sdwbar += sex_wbar*sex_wbar;            /* sum of sex_wbar^2 */
          asex_sdwbar += asex_wbar*asex_wbar;        /* sum of asex_wbar^2 */
          sex_gmwbar += log(sex_wbar);             /* sum of log(sex_wbar) */
          asex_gmwbar += log(asex_wbar);          /* sum of log(asex_wbar) */
          }

     sex_amwbar = sex_amwbar/timesup;                   /* arithmetic mean */
     asex_amwbar = asex_amwbar/timesup;
     sex_sdwbar = sqrt( (sex_sdwbar                  /* standard deviation */
                - timesup*(sex_amwbar)*(sex_amwbar))/(timesup - 1) );
     asex_sdwbar = sqrt( (asex_sdwbar
                - timesup*(asex_amwbar)*(asex_amwbar))/(timesup - 1) );
     sex_gmwbar = exp(sex_gmwbar/timesup);               /* geometric mean */
     asex_gmwbar = exp(asex_gmwbar/timesup);
     ratio_amwbar = (sex_amwbar)/(asex_amwbar);             /* comparisons */
     ratio_sdwbar = (sex_sdwbar)/(asex_sdwbar);
     ratio_gmwbar = (sex_gmwbar)/(asex_gmwbar);
     }

void init_ensemble()
     {
     double u, u2, u3, u4, v, v2, v3, v4, uv3, u2v2, u3v;
     srand48( (double) SEED);              /* seed random number generator */
     sercorr = 1.0 - p12 - p21;            /* serial correlation of states */
     probbad = p12/(p12 + p21);  /* for bad generation in stationary dist. */
     wgood = pow( (double) selection, probbad/(1.0 - probbad) );
     wbad = 1.0/((double) selection);
     u = mutrate; v = 1.0 - mutrate;
     u2 = u*u; v2 = v*v;
     u3 = u2*u; v3 = v2*v;
     u4 = u3*u; v4 = v3*v;
     uv3= u*v3; u2v2 = u2*v2; u3v = u3*v;

     mutation[0][0] = v4;          mutation[0][1] = uv3;
     mutation[0][2] = uv3;         mutation[0][3] = u2v2;
     mutation[0][4] = u2v2;        mutation[0][5] = u2v2;
     mutation[0][6] = u3v;         mutation[0][7] = u2v2;
     mutation[0][8] = u3v;         mutation[0][9] = u4;

     mutation[1][0] = 2*uv3;       mutation[1][1] = u2v2 + v4;
     mutation[1][2] = 2*u2v2;      mutation[1][3] = u3v + uv3;
     mutation[1][4] = 2*uv3;       mutation[1][5] = u3v + uv3;
     mutation[1][6] = 2*u2v2;      mutation[1][7] = 2*u3v;
     mutation[1][8] = u4 + u2v2;   mutation[1][9] = 2*u3v;

     mutation[2][0] = 2*uv3;       mutation[2][1] = 2*u2v2;
     mutation[2][2] = u2v2 + v4;   mutation[2][3] = u3v + uv3;
     mutation[2][4] = 2*u3v;       mutation[2][5] = u3v + uv3;
     mutation[2][6] = u4 + u2v2;   mutation[2][7] = 2*uv3;
     mutation[2][8] = 2*u2v2;      mutation[2][9] = 2*u3v;

     mutation[3][0] = 2*u2v2;      mutation[3][1] = u3v + uv3;
     mutation[3][2] = u3v + uv3;   mutation[3][3] = 0.5*(u4 + 2*u2v2 + v4);
     mutation[3][4] = 2*u2v2;      mutation[3][5] = 0.5*(u4 + 2*u2v2 + v4);
     mutation[3][6] = u3v + uv3;   mutation[3][7] = 2*u2v2;
     mutation[3][8] = u3v + uv3;   mutation[3][9] = 2*u2v2;

     mutation[4][0] = u2v2;        mutation[4][1] = uv3;
     mutation[4][2] = u3v;         mutation[4][3] = u2v2;
     mutation[4][4] = v4;          mutation[4][5] = u2v2;
     mutation[4][6] = uv3;         mutation[4][7] = u4;
     mutation[4][8] = u3v;         mutation[4][9] = u2v2;

     mutation[5][0] = 2*u2v2;      mutation[5][1] = u3v + uv3;
     mutation[5][2] = u3v + uv3;   mutation[5][3] = 0.5*(u4 + 2*u2v2 + v4);
     mutation[5][4] = 2*u2v2;      mutation[5][5] = 0.5*(u4 + 2*u2v2 + v4);
     mutation[5][6] = u3v + uv3;   mutation[5][7] = 2*u2v2;
     mutation[5][8] = u3v + uv3;   mutation[5][9] = 2*u2v2;

     mutation[6][0] = 2*u3v;       mutation[6][1] = 2*u2v2;
     mutation[6][2] = u4 + u2v2;   mutation[6][3] = u3v + uv3;
     mutation[6][4] = 2*uv3;       mutation[6][5] = u3v + uv3;
     mutation[6][6] = u2v2 + v4;   mutation[6][7] = 2*u3v;
     mutation[6][8] = 2*u2v2;      mutation[6][9] = 2*uv3;

     mutation[7][0] = u2v2;        mutation[7][1] = u3v;
     mutation[7][2] = uv3;         mutation[7][3] = u2v2;
     mutation[7][4] = u4;          mutation[7][5] = u2v2;
     mutation[7][6] = u3v;         mutation[7][7] = v4;
     mutation[7][8] = uv3;         mutation[7][9] = u2v2;

     mutation[8][0] = 2*u3v;       mutation[8][1] = u4 + u2v2;
     mutation[8][2] = 2*u2v2;      mutation[8][3] = u3v + uv3;
     mutation[8][4] = 2*u3v;       mutation[8][5] = u3v + uv3;
     mutation[8][6] = 2*u2v2;      mutation[8][7] = 2*uv3;
     mutation[8][8] = u2v2 + v4;   mutation[8][9] = 2*uv3;

     mutation[9][0] = u4;          mutation[9][1] = u3v;
     mutation[9][2] = u3v;         mutation[9][3] = u2v2;
     mutation[9][4] = u2v2;        mutation[9][5] = u2v2;
     mutation[9][6] = uv3;         mutation[9][7] = u2v2;
     mutation[9][8] = uv3;         mutation[9][9] = v4;

     examratio = seamratio = /* statistics for sexual to asexual ratios of */
     exsdratio = sesdratio =   /* arithmetic mean wbar, standard deviation */
     exgmratio = segmratio = 0.0;      /* of wbar, and geometric mean wbar */
     }

void updatestats(p_mean,p_standerr,p_x)  /* update intermediate statistics */
     double *p_mean, *p_standerr, *p_x; /* *p_mean used to store sum(*p_x) */
     {                          /* *p_standerr used to store sum((*p_x)^2) */
     *p_mean += *p_x;
     *p_standerr += (*p_x)*(*p_x);
     }

void finalizestats(p_mean,p_standerr,p_n)           /* finalize statistics */
     double *p_mean, *p_standerr;  /* *p_mean is estimated arithmetic mean */
     int *p_n;                /* *p_standerr is standard error of estimate */
     {          /* *p_n is sample size (ensemble size in this application) */
     *p_mean = *p_mean/(*p_n);                      /* mean, then variance */
     *p_standerr = (*p_standerr - (*p_n)*(*p_mean)*(*p_mean))/(*p_n - 1);
     *p_standerr = sqrt(*p_standerr);                /* standard deviation */
     *p_standerr = *p_standerr/sqrt((double) *p_n);      /* standard error */
     }

void put_rheadings()
     {
     printf("  Sexual Population");
     printf("  Asexual Population");
     printf("  Sex/Asex Ratio\n");
     printf("  N   AM   SD   GM");
     printf("    N   AM   SD   GM");
     printf("    AM   SD   GM\n");
     }

void put_rresults()
     {
     printf(" %3.1lf %4.2lf %4.2lf %4.2lf",
            sex_size,
            sex_amwbar,
            sex_sdwbar,
            sex_gmwbar);
     printf("  %3.1lf %4.2lf %4.2lf %4.2lf",
            asex_size,
            asex_amwbar,
            asex_sdwbar,
            asex_gmwbar);
     printf("  %4.2lf %4.2lf %4.2lf\n",
            ratio_amwbar,
            ratio_sdwbar,
            ratio_gmwbar);
     }

void put_eheadings()
     {
     printf("Select Rcmb Pbad SCor  Mutate  TimeUp Trial");
     printf(" AMwbar +/-  SDwbar +/-  GMwbar +/-\n");
     }

void put_eresults()
     {
     printf("*%5d %4.2lf %5.3lf %4.2lf %5.1le %6d %5d",
            selection,
            r,
            probbad,
            sercorr,
            mutrate,
            timesup,
            trials);
     printf(" %5.3lf %5.3lf %5.3lf %5.3lf %5.3lf %5.3lf\n\n",
            examratio,
            seamratio,
            exsdratio,
            sesdratio,
            exgmratio,
            segmratio);
     }

main()              /* main() orchestrates the runs comprising an ensemble */
{
int run;

scanf("%d",&selection);              /* read selection strength from stdin */
while (selection > 0)
     {
     scanf("%lf%lf%lf%lf%d%d",&r,&p12,&p21,&mutrate,&timesup,&trials);
     init_ensemble();                     /* initialize for a new ensemble */
     put_rheadings();
     for (run = 0; run < trials; ++run)
          {
          init_run();
          run_sex_asex();

          updatestats(&examratio,&seamratio,&ratio_amwbar);
          updatestats(&exsdratio,&sesdratio,&ratio_sdwbar);
          updatestats(&exgmratio,&segmratio,&ratio_gmwbar);
          put_rresults();
          }
     finalizestats(&examratio,&seamratio,&trials);
     finalizestats(&exsdratio,&sesdratio,&trials);
     finalizestats(&exgmratio,&segmratio,&trials);
     put_eheadings();
     put_eresults();
     scanf("%d",&selection); /* read another selection strength from stdin */
     }
}