This looks like a little complicated, but it is a very straight-forward extension of the previous genetic model.
#include <stdio.h>
#include <stdlib.h>
#define ALT 2 /* represent Altruist allele */
#define SELFISH 1 /* represent Selfish allele */
#define EMP 0 /* represent an empty space */
#define PRNT_INTVL 1 /* the interval to print out the state of population */
struct individual {
int allele;
int phenotype;
double fitness;
};
/* global parameters */
#define NUM_REPS 100
#define MAX_GEN 10000
#define POP_SIZE 200
#define INIT_FREQ 0.5 /* initial frequency of Altruist */
int gNeighborSize = 100, gDispersal=100;
double gBenefit = 0.6, gCost = 0.2;
double gColonizationDifficulty = 0.95;
double gTragedyFreq = 0;
double gDamageFraction = 0.1; /* fraction of popSize, destroyed */
int gDamageSize; /* will become popSize * gDamageFraction */
void GetPara(int *popSizePtr, double *freqPtr);
void InitPopHap(struct individual *popPtr, int popSize, double initFreq);
void CreateNextGen(struct individual *curPop, struct individual *nextPop,
int popSize, int model);
void KillEm(struct individual *pop, int popSize);
void DeterminePhenotypeHap(struct individual *popPtr, int popSize);
void AssignFitness(struct individual *popPtr, int popSize);
double NeighborAltruistFreq(struct individual *popPtr, int popSize,
int indexOfTargetIndividual, int neighborSize);
int WrappedIndex(int popSize, int index);
int RandInt(int max);
int MyRound(double x);
int CntAllele(struct individual *pop, int popSize, int allele,
int *result, int resultSize);
void PrintState(FILE *fp, struct individual *pop, int popSize);
int main(void) {
int gen, rep, popSize, numGen = MAX_GEN;
struct individual *currentPtr, *nextPtr, *popB, *popR;
int alleleCnt[3], winner[4] = {0};
double initFreq;
FILE *outFP;
srand48(time(NULL));
/* GetPara(&popSize, &initFreq); /* get the values from the user */
popSize = POP_SIZE;
initFreq = INIT_FREQ;
gDamageSize = MyRound(popSize * gDamageFraction);
/* allocate memory for the population arrays */
popB=calloc(popSize, sizeof(struct individual));
popR=calloc(popSize, sizeof(struct individual));
if (popB == NULL || popR == NULL) {
fprintf(stderr, "No Memory to allocate population arrays\n");
exit(EXIT_FAILURE);
}
/* file output of spatial structure */
outFP = fopen("altruism.spatial.out", "w");
if(outFP == NULL) {
fprintf(stderr,"Can't open the outputfile\n");
exit(EXIT_FAILURE);
}
/* start the main loop */
for (rep = 0; rep < NUM_REPS; rep++) {
InitPopHap(popB, popSize, initFreq); /* randomly assign alleles to popB */
for (gen = 0; gen < numGen; gen++) {
if (gen % 2 == 0) {
currentPtr = popB;
nextPtr = popR;
} else {
currentPtr = popR;
nextPtr = popB;
}
CreateNextGen(currentPtr, nextPtr, popSize, 1);
/* the last arg: 0=panmictic, 1=viscous */
CntAllele(nextPtr, popSize, ALT, alleleCnt, 3);
/* print changes in distribution of alleles for the first sim */
if (rep == 0 && gen % PRNT_INTVL == 0) {
PrintState(outFP, nextPtr, popSize);
}
/* Check if extinction or fixation occured */
if (alleleCnt[EMP] == popSize) { /* pop Extinct */
winner[EMP]++;
break;
} else if (alleleCnt[ALT] == 0) { /* allele selfish 0 wins */
winner[SELFISH]++;
break;
} else if (alleleCnt[SELFISH] == 0) { /* allele altruist 1 wins */
winner[ALT]++;
break;
}
}
if (gen == numGen) { /* the two alleles are segregating */
winner[3]++;
}
}
printf("Winning (extinction:selfish:alt:tie) = %d %d %d %d\n",
winner[0], winner[1], winner[2], winner[3]);
close(outFP);
free(popB);
free(popR);
return 0;
}
/* Initialize the parameters, Not complete */
void GetPara(int *popSizePtr, double *freqPtr) {
fprintf(stderr, "Population Size?: ");
scanf("%d", popSizePtr);
fprintf(stderr, "Init freq. of allele 1?: ");
scanf("%d", freqPtr);
return;
}
/*
* Assume that there are two alleles: 0 and 1
* Initialize the population, so that the allele frequency of allele 1
* is initFreq.
*/
void InitPopHap(struct individual *popPtr, int popSize, double initFreq) {
int i, numAlleles, cntr, index;
numAlleles = MyRound(popSize * initFreq);
for (i = 0; i < popSize; i++) {
popPtr[i].allele = SELFISH;
}
cntr = 0;
while (cntr < numAlleles) {
index = RandInt(popSize - 1);
if (popPtr[index].allele == SELFISH) {
popPtr[index].allele = ALT;
cntr++;
}
}
DeterminePhenotypeHap(popPtr, popSize);
AssignFitness(popPtr, popSize);
return;
}
/*
* Input:
* curPop: array containing the current population
* nextPop: array which will be created from curPop
* popSize: size of the population
* model: 0 = panmictic, 1 = viscous
*
* Side effects:
* nextPop will be created
*/
void CreateNextGen(struct individual *curPop, struct individual *nextPop,
int popSize, int model) {
int i, parentIndex;
struct individual ind;
for (i = 0; i < popSize; i ++) {
for ( ; ; ) { /* find a parent */
if (model == 0) { /* panmictic */
parentIndex = RandInt(popSize - 1);
} else { /* viscous */
parentIndex = i + RandInt(gDispersal * 2) - gDispersal;
parentIndex = WrappedIndex(popSize, parentIndex);
}
ind = curPop[parentIndex];
if (drand48() < ind.fitness) {
break;
}
}
nextPop[i] = curPop[parentIndex];
}
if (drand48() < gTragedyFreq) { /* destruction */
KillEm(nextPop, popSize);
}
DeterminePhenotypeHap(nextPop, popSize);
AssignFitness(nextPop, popSize);
return;
}
/* population disturbance */
void KillEm(struct individual *pop, int popSize) {
int i, site, index, damage;
site = RandInt(popSize - 1);
for (i = site; i < site + gDamageSize; i++){
index = WrappedIndex(popSize, i);
pop[index].allele = EMP; /* became empty */
}
return;
}
/* in the case of haploid, phenotype is exactly same as genotype */
void DeterminePhenotypeHap(struct individual *popPtr, int popSize) {
int i;
for (i = 0; i < popSize; i++) {
popPtr[i].phenotype = popPtr[i].allele;
}
return;
}
/*
* We assumed that the beginning and end of the array is connected.
* So index of -1 means the last element (popSize - 1). If the index
* is larger than popSize, e.g. popSize + 2, this correspond to the
* third individual (index of 2). It returns this "wrapped index".
*/
int WrappedIndex(int popSize, int index) {
return ((popSize + index) % popSize);
}
/*
* Input:
* popPtr: a pointer to the population array
* popSize: size of the population
*
* Side effect:
* Go through each individual in the population, calculate its fitness,
* And the value is assigned to popPtr[i].fitness.
*
* Fitness of an individual
* Altruist: 1 + b * p - c
* Selfish: 1 + b * p
* empty: d (should be less than 1 for the population to remain)
*
* Altruist pays a cost c to help others in the neighbor. The benefit
* to the neighbor group is proportional to the frequency (p) of altruist in
* the group. If everyone co-operate, the fitness of an altruist is 1 + b - c,
* so it can be higher than selfish (if benefit (b) > cost (c))
*
* The maximum individual fitness is 1 + b. So I'm actually making the
* fitness relative to 1+b. In this way fitness is between 0 and 1.
*
* The size of neighbor is determined by a global parameter gNeighborSize.
*/
void AssignFitness(struct individual *popPtr, int popSize) {
int i;
double altNeighborFreq;
for (i = 0; i < popSize; i++) {
if (popPtr[i].phenotype == EMP) { /* empty site */
popPtr[i].fitness = gColonizationDifficulty/(1+gBenefit);
} else { /* non empty site */
altNeighborFreq = NeighborAltruistFreq(popPtr, popSize, i,gNeighborSize);
if(popPtr[i].phenotype == ALT) { /* altruist fitness */
popPtr[i].fitness = (1 + altNeighborFreq*gBenefit -gCost)/(1+gBenefit);
} else if (popPtr[i].phenotype == SELFISH) { /* selfish fitness */
popPtr[i].fitness = (1 + altNeighborFreq * gBenefit) / (1+gBenefit);
}
}
}
return;
}
/*
* For the individual with indexOfTargetIndividual in a population of
* popPtr, it caclculate the freq. of altruist in the neighbors. If
* neighborSize is 1, the neighbor consists of 3 individuals (target,
* right and left ind.). With neighborSize of 2, the neighbor cosists
* of 5 individuals.
*/
double NeighborAltruistFreq(struct individual *popPtr, int popSize,
int indexOfTargetIndividual, int neighborSize){
int i, index, minIndex, maxIndex, numNeighbors, cntr;
/* set the lower and upper bound */
minIndex = indexOfTargetIndividual - neighborSize;
maxIndex = indexOfTargetIndividual + neighborSize;
/* count the individual with altruistic phenotype (ALT) */
cntr = 0;
for (i = minIndex; i <= maxIndex; i++) {
index = WrappedIndex(popSize, i);
if (popPtr[index].phenotype == ALT) {
cntr++;
}
}
numNeighbors = neighborSize * 2 + 1;
return (double) cntr / numNeighbors;
}
/*
* Returns a uniform random integer between 0 and max (ends inclusive)
*/
int RandInt(int max) {
return (int) (drand48() * (max+1));
}
/*
* Round the floating points to the nearest integer
* For example, 9.5 <= x < 10.5 become 10
* This function can cause integer over-flow when x is large
*/
int MyRound(double x) {
int intPart;
if (x > 0) {
return (int) (x + 0.5);
} else {
return (int) (x - 0.5);
}
}
/*
* Count the number of alleles of type allele in the population
*/
int CntAllele(struct individual *pop, int popSize, int allele,
int *result, int resultSize) {
int i, count=0, alleleType;
/* initialize */
for (i = 0; i < resultSize; i++) {
result[i] = 0;
}
/* count */
for (i = 0 ; i < popSize; i++) {
alleleType = pop[i].allele;
result[alleleType]++;
}
return count;
}
void PrintState(FILE *fp, struct individual *pop, int popSize) {
int i;
fprintf(fp, "%d", pop[0]);
for(i = 1; i < popSize; i++) {
fprintf(fp, ",%d", pop[i]);
}
fprintf(fp, "\n");
}