#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#define NUM_REPS 50000
const gsl_rng *gBaseRand;
double gMeanOvuleNum, gStddevOvuleNum, gSelfRate, gMeanPolVisit,
gInitFreq4x;
int gPopSize;
double NonNegNormal(double mean, double stddev);
void SetupGlobalVars(void);
unsigned int MakeNextGenSeeds (unsigned int numTargetPloidy,
unsigned int numOtherPloidy);
int CheckStatus(int gen, unsigned int num2x, unsigned int num4x);
int main(void) {
int i, rc, rep, gen, randSeed, maxGen = 10000;
unsigned int num2x, num4x, totalN;
unsigned int seeds2x, seeds4x;
double freq2x, dummyPara;
int results[4] = {0};
/* specifying to use Mersenne twister MT-19937 as the uniform PRNG */
gBaseRand = gsl_rng_alloc(gsl_rng_mt19937);
srand(time(NULL)); /* initialization for rand() */
randSeed = rand(); /* returns a non-negative integer */
gsl_rng_set (gBaseRand, randSeed); /* seed the PRNG */
SetupGlobalVars();
/* print header */
printf("Changing MeanPolVisit\n");
printf("paraVal\t4xWin\t2xWin\tPopExtinct\tCoexist\n");
for (gMeanPolVisit = 1; gMeanPolVisit < 20; gMeanPolVisit +=2) {
/* for (gMeanOvuleNum = 1; gMeanOvuleNum < 100; gMeanOvuleNum +=10) { */
/* for (gSelfRate = 0; gSelfRate <= 1; gSelfRate +=0.2) { */
/* initialize the results counter */
dummyPara = gMeanPolVisit;
for(i = 0; i < 4; i++) {
results[i] = 0;
}
for (rep = 0; rep < NUM_REPS; rep++) {
/* beginning of a single simulation */
num2x = (int) (gPopSize * (1-gInitFreq4x));
num4x = (int) (gPopSize * gInitFreq4x);
for (gen = 0; gen < maxGen; gen++) {
seeds2x = MakeNextGenSeeds(num2x, num4x);
seeds4x = MakeNextGenSeeds(num4x, num2x);
if (seeds2x + seeds4x > 0) {
freq2x = (double) seeds2x / (seeds2x + seeds4x);
num2x = gsl_ran_binomial(gBaseRand, freq2x, gPopSize);
num4x = gPopSize - num2x;
} else {
num2x = num4x = 0;
}
rc = CheckStatus(gen, num2x, num4x);
if (rc > 0) { /* one ploidy went extinct */
results[rc - 1]++;
break;
}
}
if (gen == maxGen) { /* coexistence */
results[3]++;
}
}
printf("%lf\t%d\t%d\t%d\t%d\n", dummyPara,
results[0], results[1], results[2], results[3]);
}
gsl_rng_free(gBaseRand);
return 0;
}
/*
* This function returns a random deviate sampled from truncated
* normal distribution (i.e., the value returned is non negative).
*/
double NonNegNormal(double mean, double stddev) {
double normalDeviate = -1.0;
while (normalDeviate < 0) {
normalDeviate = gsl_ran_gaussian (gBaseRand, stddev) + mean;
}
return normalDeviate;
}
/*
* set up following variables: gMeanOvuleNum, gStddevOvuleNum,
* gSelfRate, gMeanPolVisit;
*/
void SetupGlobalVars(void) {
gMeanOvuleNum = 10;
gStddevOvuleNum = gMeanOvuleNum/2;
gSelfRate = 0.8;
gMeanPolVisit = 3;
gInitFreq4x = 0.2;
gPopSize = 20;
/* use the following for interactive setup
fprintf(stderr, "Population Size: ");
scanf("%d", &gPopSize);
fprintf(stderr, "Init freq. of tetraploids: ");
scanf("%lf", &gInitFreq4x);
fprintf(stderr, "Mean number of ovules: ");
scanf("%lf", &gMeanOvuleNum);
fprintf(stderr, "Std. dev. of ovule numbers: ");
scanf("%lf", &gStddevOvuleNum);
fprintf(stderr, "Selfing rate [0 - 1]: ");
scanf("%lf", &gSelfRate);
fprintf(stderr, "Mean number of pollinator visit per flower: ");
scanf("%lf", &gMeanPolVisit);
*/
printf ("Parameters: meanOvNum = %lf, stddevOvNum = %lf, selfRate = %lf, "
"meanPolVisit = %lf, initNum2x = %d, initNum4x = %d\n",
gMeanOvuleNum, gStddevOvuleNum, gSelfRate,
gMeanPolVisit,
(int) (gPopSize * (1-gInitFreq4x)), (int) (gPopSize * gInitFreq4x));
return;
}
/*
* Take the state of current population as the argument.
*
* Return the number of seeds of the target ploidy for the next generation.
*
* For example, if the 1st argument is the number of diploid individuals and
* the 2nd argument is the number of tetraploid individuals, it goes through
* the mating process, and returns the number of diploid seeds in the
* next generation. If the 2 arguments are switched (i.e. numTargetPloidy
* is the number of tetraploid individuals), it will return the
* number of tetraploid seeds in the next generation.
*/
unsigned int MakeNextGenSeeds (unsigned int numTargetPloidy,
unsigned int numOtherPloidy) {
int i;
unsigned int numOv, numSelfedOv, numOCOv, numPollinatorVisit, numPollen;
unsigned int goodOCSeeds, sum = 0;
double freqThisPloidy, freqThisPloidyInPollen;
/* frequency of this ploidy in the adults */
freqThisPloidy = (double) numTargetPloidy/(numTargetPloidy + numOtherPloidy);
/* go through each individual and let it produce seeds */
for (i = 0; i < numTargetPloidy; i++) {
/* number of ovules which this individual produces */
numOv = (int) NonNegNormal(gMeanOvuleNum, gStddevOvuleNum);
/* some of these ovules are used for selfing */
numSelfedOv = gsl_ran_binomial(gBaseRand, gSelfRate, numOv);
/* remainders are available for outcrossing */
numOCOv = numOv - numSelfedOv;
/* number of pollinators visitied this flower */
numPollinatorVisit = gsl_ran_poisson(gBaseRand, gMeanPolVisit);
/* number of outcrossed seeds which received the correct ploidy pollen */
if (numPollinatorVisit > 0) {
/* Freq. of pollen with the same ploidy as the mom */
freqThisPloidyInPollen =
(double) gsl_ran_binomial(gBaseRand,freqThisPloidy,numPollinatorVisit)/
numPollinatorVisit;
goodOCSeeds =
gsl_ran_binomial(gBaseRand, freqThisPloidyInPollen, numOCOv);
} else { /* no pollinator visit, so no outcrossed seeds */
goodOCSeeds = 0;
}
sum += numSelfedOv + goodOCSeeds; /* tracking the total viable seeds */
}
return sum;
}
#define PRINT_INTVL 10
/* This function returns
* 0: 2 ploidy exists
* 1: diploid went extinct
* 2: tetraploid went extinct
* 3: population went extinct
*/
int CheckStatus(int gen, unsigned int num2x, unsigned int num4x) {
/* if (gen % PRINT_INTVL == 0) {
printf ("%d\t%u\t%u\n", gen, num2x, num4x);
}
*/
if (num2x == 0 && num4x == 0){
return 3;
}
if (num2x == 0) {
return 1;
}
if (num4x == 0) {
return 2;
}
return 0;
}