// BottleSim source code // Chih-Horng Kuo // preprocessor directive #include #include #include #include #include #include #include #include #include #include #include // For standard Template Library (STL) using namespace std; // constants const char programName[21] = "BottleSim"; // program name const char versionNumber[15] = "2.6.1(Linux)"; // version number const char versionDate[15] = "01/09/2006"; // version date const char author[61] = "Chih-Horng Kuo & Fred Janzen"; // author const char email[21] = "chkuo@lifedev.org"; // author email const char web[61] = "http://chkuo.name/"; // author web const int fileNameLength = 21; // filename length const int locusNameLength = 11; // locus name length const int idLength = 11; // length of character array that holds individual id const char inputMissingValue = '?'; // default value for input missing value const short maxAllele = 62; // maximum number of alleles per locus // possible allele list for converting genotype input file const char arrayAlleleId[maxAllele] = {'0', '1', '2', '3', '4', '5', '6', '7', '8', '9', 'A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'J', 'K', 'L', 'M', 'N', 'O', 'P', 'Q', 'R', 'S', 'T', 'U', 'V', 'W', 'X', 'Y', 'Z', 'a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'i', 'j', 'k', 'l', 'm', 'n', 'o', 'p', 'q', 'r', 's', 't', 'u', 'v', 'w', 'x', 'y', 'z'}; // typedef typedef char fileName[fileNameLength]; // new data type "fileName" to hold file name typedef char locusName[locusNameLength]; // new data type "locusName" to hold locus name typedef char individualId[idLength]; // new data type "individualId" to hold individual id // function prototypes // infomation display void welcome(void); void memoryError(void); // get settings int moduleSelection(void); void getCPSSimSetting(int °reeOverlap, int &reproSys, int &longevity, int &ageMaturation, int &nPB, int &nB, int &nFemale, int &nMale, int &nYear, int &nIteration, int &boolOutputRawGenotype); void getVPSSimSetting(int °reeOverlap, int &reproSys, int &longevity, int &ageMaturation, int &nIteration, int &boolOutputRawGenotype); int getDegreeOverlap(void); int getReproSys(void); int getLongevity(void); int getAgeMaturation(int longevity); int getNPB(void); int getNB(void); int getNFemale(int nB); int getNYear(void); int getNIteration(void); int getBoolOutputRawGenotype(void); int getNMaxPopSize(int *nPopSize, int nYear); // get random number int getRandomNumber(int min, int max); // get allele count void getDiSLAlleleCount(int *arrayAlleleCount, double *inputAlleleFreq, int nAllele, int nPopSize, int allelePerLocus); // get new genotype void getNewDiSLGenotype(int **arrayParent, int *arrayNewGenotype, int reproSys, int nAdultTotal, int nFemale, int parent1, int parent2); void getNewDiMLGenotype(int ***arrayParent, int **arrayNewGenotype, int reproSys, int nAdultTotal, int nFemale, int parent1, int parent2, int nLoci); // read input file int readInputFileNAllele(ifstream &alleleFreqInput); void readInputFileAlleleFreq(ifstream &alleleFreqInput, double *inputAlleleFreq, int nAllele); int readInputFileNYear(ifstream &popSizeInput); void readInputFilePopSize(ifstream &popSizeInput, int *nPopSize, int *nFemale, int *nMale, int nYear, int reproSys); int readInputFileNSize(ifstream &genotypeInputFile); int readInputFileNLoci(ifstream &genotypeInputFile); void readInputFileLocusName(ifstream &genotypeInputFile, locusName *arrayInputLocus, int nLoci); void readInputFileGenotype(ifstream &genotypeInputFile, individualId *arrayInputId, char ***arrayInput, int nSize, int nLoci, int allelePerLocus); // display void displayAlleleFreqInput(double *inputAlleleFreq, int nAllele); void displayVPopSizeInput(int *nPopSize, int *nFemale, int reproSys, int nYear); void displayGenotypeInput(locusName *arrayInputLocus, individualId *arrayInputId, char ***arrayInput, int nSize, int nLoci, int allelePerLocus); // verify input void checkAlleleFreqInput(double *inputAlleleFreq, int nAllele); void checkVPopSsizeInput(int *nPopSize, int *nFemale, int reproSys, int nYear); // write setting void writeOutputHeader(ofstream &outputFile); void writeSLRawGenotypeHeader(ofstream &rawGenotypeOutputFile); void writeMLRawGenotypeHeader(ofstream &rawGenotypeOutputFile, locusName *arrayInputLocus, int nLoci); void writeMainSetting(ofstream &outputFile, int module, fileName inputFileName, fileName outputFileName); void writeVPopSizeInput(ofstream &outputFile, int *nPopSize, int *nFemale, int reproSys, int nYear); void writeAlleleFreq(ofstream &outputFile, int nLoci, locusName *arrayInputLocus, float **arrayAlleleFreq); void writeCPSSimSetting (ofstream &outputFile, int degreeOverlap, int reproSys, int longevity, int ageMaturation, int nPB, int nB, int nFemale, int nMale, int nYear, int nIteration); void writeVPSSimSetting (ofstream &outputFile, int degreeOverlap, int reproSys, int longevity, int ageMaturation, int nYear, int nIteration); void writeDegreeOverlap(ofstream &outputFile, int degreeOverlap); void writeReprosys(ofstream &outputFile, int reproSys); void writeLongevity(ofstream &outputFile, int longevity); void writeAgeMaturation(ofstream &outputFile, int ageMaturation); void writeNPB(ofstream &outputFile, int nPB); void writeNB(ofstream &outputFile, int nB); void writeSexRatio(ofstream &outputFile, int nFemale, int nMale); void writeNYear(ofstream &outputFile, int nYear); void writeNIteration(ofstream &outputFile, int nIteration); // write output void writeDiSLSummary(ofstream &outputFile, int nYear, int nIteration, int **arrayRawOA, float **arrayRawEA, float **arrayRawHo, float **arrayRawHe, float **arrayRawF); void writeDiMLSummary(ofstream &outputFile, int nLoci, locusName *arrayInputLocus, int nYear, int nIteration, int ***arrayRawOA, float ***arrayRawEA, float ***arrayRawHo, float ***arrayRawHe, float ***arrayRawF); void writeFixDiploidSingleLocus(ofstream &outputFile, int nIteration, int fixCase, float fixProb); void writeFixDiploidMultiLocus(ofstream &outputFile, int nLoci, locusName *arrayInputLocus, int nIteration, int *fixCase, float *fixProb); void writeSingleLocus(ofstream &outputFile, int nYear, float *arrayMean, float *arrayStdD, float *arrayStdE, float *arrayPercent); void writeMultiLocus(ofstream &outputFile, int nLoci, locusName *arrayInputLocus, int nYear, float **arrayMean, float *arrayAv, float *arrayStdD, float *arrayStdE, float *arrayPercent); void writeDiSLRawGenotype(ofstream &rawGenotypeOutputFile, int **arrayGenotype, int allelePerLocus, int nPopSize); void writeDiMLRawGenotype(ofstream &rawGenotypeOutputFile, int ***arrayGenotype, int allelePerLocus, int nPopSize, int nLoci); // calculation void statGenotypeInput(char ***arrayInput, int **arrayAlleleCount, int *alleleTotalCount, float **arrayAlleleFreq, int nSize, int nLoci, int allelePerLocus); void statFixDiploidSingleLocus(int **arrayRawOA, int nYear, int nIteration, int &fixCase, float &fixProb); void statFixDiploidMultiLocus(int ***arrayRawOA, int nLoci, int nYear, int nIteration, int *fixCase, float *fixProb); void statDiploidSingleLocus(int **arrayGenotype, int allelePerLocus, int nPopSize, int nAllele, int &nOA, float &nEA, float &nHo, float &nHe, float &nF); void statDiploidMultiLocus(int ***arrayGenotype, int allelePerLocus, int nPopSize, int nLoci, int *arrayOA, float *arrayEA, float *arrayHo, float *arrayHe, float *arrayF); void outputDiploidSingleLocus(int **arrayRawdata, int nYear, int nIteration, float *arrayMean, float *arrayStdD, float *arrayStdE, float *arrayPercent); void outputDiploidSingleLocus(float **arrayRawdata, int nYear, int nIteration, float *arrayMean, float *arrayStdD, float *arrayStdE, float *arrayPercent); void outputDiploidMultiLocus(int ***arrayRawdata, int nLoci, int nYear, int nIteration, float **arrayMean, float *arrayAv, float *arrayStdD, float *arrayStdE, float *arrayPercent); void outputDiploidMultiLocus(float ***arrayRawdata, int nLoci, int nYear, int nIteration, float **arrayMean, float *arrayAv, float *arrayStdD, float *arrayStdE, float *arrayPercent); // file operation void openInput(ifstream &inputFile, fileName &inputFileName); void openPopSizeInput(ifstream &popSizeInput, fileName &popSizeInputName); void openOutput(ofstream &outputFile, fileName &outputFileName); void openRawGenotypeOutput(ofstream &rawGenotypeOutputFile, fileName &rawGenotypeOutputName); void closeInputFile(ifstream &inputFile); void closeOutputFile(ofstream &outputFile); // simulation module void simDiploidSingleLocusCPS(void); void simDiploidSingleLocusVPS(void); void simDiploidMultiLocusCPS(void); void simDiploidMultiLocusVPS(void); // start of program functions // main function int main(void) { // call welcome function to display welcome message welcome(); // seed random number generator with current system time srand(time(0)); // declare variables int module = 0; //module used for the run // call moduleSelection function module = moduleSelection(); // call selected module switch (module) { case 1: // diploid singlelocus, constant population size // call simDiploidSingleLocusCPS function simDiploidSingleLocusCPS(); break; case 2: // diploid singlelocus, variable population size // call simDiploidSingleLocusVPS function simDiploidSingleLocusVPS(); break; case 3: // diploid multilocus, constant population size // call simDiploidMultiLocusCPS function simDiploidMultiLocusCPS(); break; case 4: // diploid multilocus, variable population size // call simDiploidMultiLocusVPS function simDiploidMultiLocusVPS(); break; default: cout << endl; cout << "*** Simulation module selection error" << endl; exit(0); } // end of switch // end, display message cout << endl; cout << "*** Program ended. ***" << endl; return 0; } // end of main function // welcome function // display welcome info void welcome(void) { // display information for users cout << "*** " << programName << " ***" << endl << endl; cout << "Version number: " << versionNumber << endl; cout << "Release date: " << versionDate << endl; cout << endl; cout << "Contact information" << endl; cout << " Author: " << author << endl; cout << " email: " << email << endl; cout << " Web: " << web << endl; cout << endl; cout << "Description:" << endl; cout << " This program is a population bottleneck simulator. " << "Please see distibution web site for detailed description, " << "input file format, program limitation, " << "and version updates." << endl; cout << endl; } // end of welcome function // memoryError function void memoryError(void) { cout << endl; cout << "*** Error allocating memory!" << endl; exit(0); } // moduleSelection function // display modules available and return user selection to main function int moduleSelection(void) { // declare variables int module = 0; // display available modules cout << endl; cout << "Simulation modules available:" << endl; cout << " (1) Diploid singlelocus, constant population size" << endl; cout << " (2) Diploid singlelocus, variable population size" << endl; cout << " (3) Diploid multilocus, constant population size" << endl; cout << " (4) Diploid multilocus, variable population size" << endl; // prompt for user selection do { cout << endl; cout << "Please choose simulation module: "; cin >> module; } while (module < 1 || module > 4); // end of do-while loop // return user selection to main cout << endl; cout << "Module selected = " << module << endl; return module; } // end of moduleSelection function // getCPSSimSetting function void getCPSSimSetting(int °reeOverlap, int &reproSys, int &longevity, int &ageMaturation, int &nPB, int &nB, int &nFemale, int &nMale, int &nYear, int &nIteration, int &boolOutputRawGenotype) { // get simulation settings // display cout << endl; cout << "Simulation settings:" << endl; // call get simulation setting functions degreeOverlap = getDegreeOverlap(); // overlapping generation setting reproSys = getReproSys(); // reproductive system option longevity = getLongevity(); // longevity the organism ageMaturation = getAgeMaturation(longevity); // age of reproduvtive maturation nPB = getNPB(); // pre-bottleneck population size nB = getNB(); // bottlneck population size // number of female/male for diecious models if (reproSys == 5 || reproSys == 6 || reproSys == 7) { nFemale = getNFemale(nB); // number of female individuals nMale = (nB - nFemale); // number of male individuals } nYear = getNYear(); // number of years to simulate nIteration = getNIteration(); // number of iterations of the run boolOutputRawGenotype = getBoolOutputRawGenotype(); // bool for output raw genotype data } // end of getCPSSimSetting function // getVPSSimSetting function void getVPSSimSetting(int °reeOverlap, int &reproSys, int &longevity, int &ageMaturation, int &nIteration, int &boolOutputRawGenotype) { // get simulation settings // display cout << endl; cout << "Simulation settings:" << endl; // call get simulation setting functions degreeOverlap = getDegreeOverlap(); // overlapping generation setting reproSys = getReproSys(); // reproductive system option longevity = getLongevity(); // longevity of the organism ageMaturation = getAgeMaturation(longevity); // age of reproduvtive maturation nIteration = getNIteration(); // number of iterations of the run boolOutputRawGenotype = getBoolOutputRawGenotype(); // bool for output raw genotype data } // end of getVPSSimSetting function // getDegreeOverlap function int getDegreeOverlap(void) { int degreeOverlap; do { cout << endl << "(0) = Minimum overlap" << endl << " (discrete generations, all individual start with age 0)" << endl << "(100) = Maximum overlap" << endl << " (all individual start with random age value)" << endl << "Please enter generation overlap setting (0-100): "; cin >> degreeOverlap; } while (degreeOverlap < 0 || degreeOverlap > 100); return degreeOverlap; } // end of getDegreeOverlap function // getReproSys function int getReproSys(void) { int reproSys; do { cout << endl << "(1) = Asexual reproduction" << endl << "(2) = Monoecy with complete selfing" << endl << "(3) = Monoecy with random mating (with selfing)" << endl << "(4) = Monoecy with random mating (without selfing)" << endl << "(5) = Dioecy with random mating" << endl << "(6) = Dioecy with single reproducing male each year" << endl << "(7) = Dioecy with single reproducing pair each year" << endl << "Please choose reproductive system of the organism: "; cin >> reproSys; } while (reproSys < 1 || reproSys > 7); return reproSys; } // end of getReproSys function // getLongevity function int getLongevity(void) { int longevity; do { cout << endl << "Please enter the longevity of the organism (year): "; cin >> longevity; } while (longevity < 1); return longevity; } // end of getLongevity function // getAgeMaturation function int getAgeMaturation(int longevity) { int ageMaturation; do { cout << endl << "Please enter the age of reproductive maturation: "; cin >> ageMaturation; } while (ageMaturation < 0 || ageMaturation > longevity); return ageMaturation; } // end of getAgeMaturation function // getNPB function int getNPB(void) { int nPB; do { cout << endl << "Please enter the population size before the bottleneck: "; cin >> nPB; } while (nPB < 1); return nPB; } // end of getNPB function // getNB function int getNB(void) { int nB; do { cout << endl << "Please enter the population size during the bottleneck: "; cin >> nB; } while (nB < 1); return nB; } // end of getNB function // getNFemale function int getNFemale(int nB) { int nFemale; do { cout << endl << "Please enter the number of females in the population during bottleneck: "; cin >> nFemale; } while (nFemale < 1 || nFemale > nB); return nFemale; } // end of getNFemale function // getNYear function int getNYear(void) { int nYear; do { cout << endl << "Please enter the number of years to simulate: "; cin >> nYear; } while (nYear < 0); return nYear; } // end of getNYear function // getNIteration function int getNIteration(void) { int nIteration; do { cout << endl << "Please enter the number of iterations you want to perform: "; cin >> nIteration; } while (nIteration < 1); return nIteration; } // end of getNIteration function // getBoolOutputRawGenotype function int getBoolOutputRawGenotype(void) { int boolOutputRawGenotype; char yesNo; do { cout << endl << "Generate genotype raw data output file?" << endl << " ***CAUTION***" << endl << " The raw data output file can be very large, " << endl << " please consider the simulation settings used" << endl << " (population size, loci number, iteration number)." << endl << endl << "Enter Y or N: "; cin >> yesNo; } while (yesNo != 'Y' && yesNo != 'y' && yesNo != 'N' && yesNo != 'n'); if (yesNo == 'Y' || yesNo == 'y') boolOutputRawGenotype = 1; else boolOutputRawGenotype = 0; return boolOutputRawGenotype; } // end of getBoolOutputRawGenotype function // getNMaxPopSize function int getNMaxPopSize(int *nPopSize, int nYear) { int nMaxPopSize; // maximum population size in pop size input int countY; nMaxPopSize = nPopSize[0]; for (countY = 0; countY < (nYear+1); countY++) { if (nMaxPopSize < nPopSize[countY]) { nMaxPopSize = nPopSize[countY]; } } return nMaxPopSize; } // end of getNMaxPopSize function // getRandomNumber function // receive minimum and maximum of the random number desired // return random number within the min/max range int getRandomNumber(int min, int max) { return ((rand() % (max + 1 - min)) + min); } // end of getRandomNumber function // getDiSLAlleleCount function void getDiSLAlleleCount(int *arrayAlleleCount, double *inputAlleleFreq, int nAllele, int nPopSize, int allelePerLocus) { int countA; // calculate allele count for each allele according to input allele freq for (countA = 0; countA < nAllele; countA++) { arrayAlleleCount[countA] = int ((nPopSize * allelePerLocus) * inputAlleleFreq[countA]); } // check allele count int sumAlleleCount = 0; for (countA = 0; countA < nAllele; countA++) { sumAlleleCount += arrayAlleleCount[countA]; } if (sumAlleleCount != (nPopSize * allelePerLocus)) { cout << endl; cout << "*** Input allele frequency error." << endl; cout << "*** Allele frequency in the input file resulted in non-integer allele count in year0 population." << endl; exit(0); } } // end of getDiSLAlleleCount function // getNewDiSLGenotype function void getNewDiSLGenotype(int **arrayParent, int *arrayNewGenotype, int reproSys, int nAdultTotal, int nFemale, int parent1, int parent2) { // check nAdultTotal/nFemale/nMale switch (reproSys) { case 1: // asexual reproduction case 2: // Monoecious - complete selfing case 3: // Monoecious - random mating with selfing if (nAdultTotal < 1) // check nAdultTotal { cout << "*** Population size setting error" << endl << "*** Number of adults < 1" << endl << "*** Unable to generate new individuals" << endl; exit(0); } break; case 4: // Monoecious - random mating without selfing if (nAdultTotal < 2) // check nAdultTotal { cout << "*** Population size setting error" << endl << "*** Number of adults < 2" << endl << "*** Unable to generate new individuals" << endl; exit(0); exit(0); } break; case 5: // Diecious - random mating case 6: // Diecious - single reproducing male case 7: // Diecious - single reproducing pair if (nAdultTotal < 2) // check nAdultTotal { cout << "*** Population size setting error" << endl << "*** Number of adults < 2" << endl << "*** Unable to generate new individuals" << endl; exit(0); exit(0); } if (nFemale < 1) // check nFemale { cout << "*** Population size setting error" << endl << "*** Number of breeding female < 1" << endl << "*** Unable to generate new individuals" << endl; exit(0); } if ((nAdultTotal - nFemale) < 1) // check nMale { cout << "*** Population size setting error" << endl << "*** Number of breeding male < 1" << endl << "*** Unable to generate new individuals" << endl; exit(0); } break; default: cout << endl; cout << "*** Reproductive system setting error" << endl; exit(0); } // end of reproSys switch // get new genotype switch (reproSys) { case 1: // asexual reproduction // get parent1 parent1 = getRandomNumber(0, (nAdultTotal - 1)); // allele1 from parent1 arrayNewGenotype[0] = arrayParent[parent1][0]; // allele2 from parent1 arrayNewGenotype[1] = arrayParent[parent1][1]; break; case 2: // Monoecious - complete selfing // get parent1 parent1 = getRandomNumber(0, (nAdultTotal - 1)); // allele1 from parent1 arrayNewGenotype[0] = arrayParent[parent1][getRandomNumber(0, 1)]; // allele2 from parent1 arrayNewGenotype[1] = arrayParent[parent1][getRandomNumber(0, 1)]; break; case 3: // Monoecious - random mating with selfing // get parent1 parent1 = getRandomNumber(0, (nAdultTotal - 1)); // get parent2 parent2 = getRandomNumber(0, (nAdultTotal - 1)); // allele from parent1 arrayNewGenotype[0] = arrayParent[parent1][getRandomNumber(0, 1)]; // allele from parent2 arrayNewGenotype[1] = arrayParent[parent2][getRandomNumber(0, 1)]; break; case 4: // Monoecious - random mating without selfing // get parent1 parent1 = getRandomNumber(0, (nAdultTotal - 1)); // get parent2 do { parent2 = getRandomNumber(0, (nAdultTotal - 1)); } while (parent2 == parent1); // allele from parent1 arrayNewGenotype[0] = arrayParent[parent1][getRandomNumber(0, 1)]; // allele from parent2 arrayNewGenotype[1] = arrayParent[parent2][getRandomNumber(0, 1)]; break; case 5: // Diecious - random mating // get parent1 and parent2 parent1 = getRandomNumber(0, (nFemale - 1)); parent2 = getRandomNumber(nFemale, (nAdultTotal - 1)); // allele from parent1 arrayNewGenotype[0] = arrayParent[parent1][getRandomNumber(0, 1)]; // allele from parent2 arrayNewGenotype[1] = arrayParent[parent2][getRandomNumber(0, 1)]; break; case 6: // Diecious - single reproducing male // get parent1 parent1 = getRandomNumber(0, (nFemale - 1)); // allele from parent1 arrayNewGenotype[0] = arrayParent[parent1][getRandomNumber(0, 1)]; // allele from parent2 arrayNewGenotype[1] = arrayParent[parent2][getRandomNumber(0, 1)]; break; case 7: // Diecious - single reproducing pair // allele from parent1 arrayNewGenotype[0] = arrayParent[parent1][getRandomNumber(0, 1)]; // allele from parent2 arrayNewGenotype[1] = arrayParent[parent2][getRandomNumber(0, 1)]; break; default: cout << endl; cout << "*** Reproductive system setting error" << endl; exit(0); } // end of reproSys switch } // end of getNewDiSLGenotype // getNewDiMLGenotype function void getNewDiMLGenotype(int ***arrayParent, int **arrayNewGenotype, int reproSys, int nAdultTotal, int nFemale, int parent1, int parent2, int nLoci) { int countL; // check nAdultTotal/nFemale/nMale switch (reproSys) { case 1: // asexual reproduction case 2: // Monoecious - complete selfing case 3: // Monoecious - random mating with selfing if (nAdultTotal < 1) // check nAdultTotal { cout << "*** Population size setting error" << endl << "*** Number of adults < 1" << endl << "*** Unable to generate new individuals" << endl; exit(0); } break; case 4: // Monoecious - random mating without selfing if (nAdultTotal < 2) // check nAdultTotal { cout << "*** Population size setting error" << endl << "*** Number of adults < 2" << endl << "*** Unable to generate new individuals" << endl; exit(0); exit(0); } break; case 5: // Diecious - random mating case 6: // Diecious - single reproducing male case 7: // Diecious - single reproducing pair if (nAdultTotal < 2) // check nAdultTotal { cout << "*** Population size setting error" << endl << "*** Number of adults < 2" << endl << "*** Unable to generate new individuals" << endl; exit(0); exit(0); } if (nFemale < 1) // check nFemale { cout << "*** Population size setting error" << endl << "*** Number of breeding female < 1" << endl << "*** Unable to generate new individuals" << endl; exit(0); } if ((nAdultTotal - nFemale) < 1) // check nMale { cout << "*** Population size setting error" << endl << "*** Number of breeding male < 1" << endl << "*** Unable to generate new individuals" << endl; exit(0); } break; default: cout << endl; cout << "*** Reproductive system setting error" << endl; exit(0); } // end of reproSys switch // get new genotype switch (reproSys) { case 1: // asexual reproduction // get parent1 parent1 = getRandomNumber(0, (nAdultTotal - 1)); for (countL = 0; countL < nLoci; countL++) { // allele1 from parent1 arrayNewGenotype[countL][0] = arrayParent[parent1][countL][0]; // allele2 from parent1 arrayNewGenotype[countL][1] = arrayParent[parent1][countL][1]; } break; case 2: // Monoecious - complete selfing // get parent1 parent1 = getRandomNumber(0, (nAdultTotal - 1)); for (countL = 0; countL < nLoci; countL++) { // allele1 from parent1 arrayNewGenotype[countL][0] = arrayParent[parent1][countL][getRandomNumber(0, 1)]; // allele2 from parent1 arrayNewGenotype[countL][1] = arrayParent[parent1][countL][getRandomNumber(0, 1)]; } break; case 3: // Monoecious - random mating with selfing // get parent1 parent1 = getRandomNumber(0, (nAdultTotal - 1)); // get parent2 parent2 = getRandomNumber(0, (nAdultTotal - 1)); for (countL = 0; countL < nLoci; countL++) { // allele from parent1 arrayNewGenotype[countL][0] = arrayParent[parent1][countL][getRandomNumber(0, 1)]; // allele from parent2 arrayNewGenotype[countL][1] = arrayParent[parent2][countL][getRandomNumber(0, 1)]; } break; case 4: // Monoecious - random mating without selfing // get parent1 and parent2 parent1 = getRandomNumber(0, (nAdultTotal - 1)); // get parent 2 do { parent2 = getRandomNumber(0, (nAdultTotal - 1)); } while (parent2 == parent1); for (countL = 0; countL < nLoci; countL++) { // allele from parent1 arrayNewGenotype[countL][0] = arrayParent[parent1][countL][getRandomNumber(0, 1)]; // allele from parent2 arrayNewGenotype[countL][1] = arrayParent[parent2][countL][getRandomNumber(0, 1)]; } break; case 5: // Diecious - random mating // get parent1 and parent2 parent1 = getRandomNumber(0, (nFemale - 1)); parent2 = getRandomNumber(nFemale, (nAdultTotal - 1)); for (countL = 0; countL < nLoci; countL++) { // allele from parent1 arrayNewGenotype[countL][0] = arrayParent[parent1][countL][getRandomNumber(0, 1)]; // allele from parent2 arrayNewGenotype[countL][1] = arrayParent[parent2][countL][getRandomNumber(0, 1)]; } break; case 6: // Diecious - single reproducing male // get parent1 parent1 = getRandomNumber(0, (nFemale - 1)); for (countL = 0; countL < nLoci; countL++) { // allele from parent1 arrayNewGenotype[countL][0] = arrayParent[parent1][countL][getRandomNumber(0, 1)]; // allele from parent2 arrayNewGenotype[countL][1] = arrayParent[parent2][countL][getRandomNumber(0, 1)]; } break; case 7: // Diecious - single reproducing pair for (countL = 0; countL < nLoci; countL++) { // allele from parent1 arrayNewGenotype[countL][0] = arrayParent[parent1][countL][getRandomNumber(0, 1)]; // allele from parent2 arrayNewGenotype[countL][1] = arrayParent[parent2][countL][getRandomNumber(0, 1)]; } break; default: cout << endl; cout << "*** Reproductive system setting error" << endl; exit(0); } // end of reproSys switch } // end of getNewDiMLGenotype // readInputFileNAllele function int readInputFileNAllele(ifstream &alleleFreqInput) { int nAllele; alleleFreqInput >> nAllele; if (nAllele < 1) { cout << endl; cout << "*** Observed number of alleles < 1" << endl; exit(0); } return nAllele; } // readInputFileAlleleFreq function void readInputFileAlleleFreq(ifstream &alleleFreqInput, double *inputAlleleFreq, int nAllele) { int countA; for (countA = 0; countA < nAllele; countA++) { if ( alleleFreqInput.eof()) break; alleleFreqInput >> inputAlleleFreq[countA]; } // display cout << "Finished reading allele frequency input file." << endl; cout << endl; } // end of readInputFileAlleleFreq function // readInputFileNYear function int readInputFileNYear(ifstream &popSizeInput) { int nYear; popSizeInput >> nYear; if (nYear < 0) { cout << endl; cout << "*** Population size input error." << endl; cout << "*** Number of years < 0" << endl; exit(0); } return nYear; } // end of readInputFileNYear function // readInputFilePopSize function void readInputFilePopSize(ifstream &popSizeInput, int *nPopSize, int *nFemale, int *nMale, int nYear, int reproSys) { int countY; // loop counting int dummyYear; // read the "year" value in the input file, do NOT actully use the number for (countY = 0; countY < (nYear+1); countY++) { if (popSizeInput.eof()) break; popSizeInput >> dummyYear; popSizeInput >> nPopSize[countY]; if (reproSys >= 5 && reproSys <= 7) { popSizeInput >> nFemale[countY]; // read nFemale nMale[countY] = (nPopSize[countY] - nFemale[countY]); // calculate nMale } } // display cout << "Finished reading population size input file." << endl; cout << endl; } // end of readInputFilePopSize function // readInputFileNSize function int readInputFileNSize(ifstream &genotypeInputFile) { int nSize; genotypeInputFile >> nSize; if (nSize < 1) { cout << endl; cout << "*** Number of individuals in the input file < 1" << endl; exit(0); } return nSize; } // end of readInputFileNSize function // readInputFileNLoci function int readInputFileNLoci(ifstream &genotypeInputFile) { int nLoci; genotypeInputFile >> nLoci; if (nLoci < 1) { cout << endl; cout << "*** Number of loci in the input file < 1" << endl; exit(0); } return nLoci; } // end of readInputFileNLoci function // readInputFileLocusName function void readInputFileLocusName(ifstream &genotypeInputFile, locusName *arrayInputLocus, int nLoci) { int countL; for ( countL = 0; countL < nLoci; countL++) { genotypeInputFile >> setw(locusNameLength) >> arrayInputLocus[countL]; if ( genotypeInputFile.eof()) break; } } // end of readInputFileLocusName function // readInputFileGenotype function void readInputFileGenotype(ifstream &genotypeInputFile, individualId *arrayInputId, char ***arrayInput, int nSize, int nLoci, int allelePerLocus) { int countN; int countL; int countA; for (countN = 0; countN < nSize; countN++) { genotypeInputFile >> setw(idLength) >> arrayInputId[countN]; if ( genotypeInputFile.eof()) break; for ( countL = 0; countL < nLoci; countL++) { for ( countA = 0; countA < allelePerLocus; countA++) { genotypeInputFile >> arrayInput[countN][countL][countA]; } } } // display cout << "Finished reading input file." << endl; cout << endl; } // end of readInputFileGenotype function // displayAlleleFreqInput function void displayAlleleFreqInput(double *inputAlleleFreq, int nAllele) { int countA; cout << "Imported data from input file: " << endl; cout << " Total number of alleles = " << nAllele << endl; for ( countA = 0; countA < nAllele; countA++) { cout << " Frequency of allele " << (countA + 1) << " = " << setprecision(4) << setiosflags(ios::fixed | ios::showpoint) << inputAlleleFreq[countA] << endl; } cout << endl; } // end of displayAlleleFreqInput function // displayVPopSizeInput function void displayVPopSizeInput(int *nPopSize, int *nFemale, int reproSys, int nYear) { int countY; cout << "Imported data from population size input file: " << endl; cout << " Number of years to simulate = " << nYear << endl; cout << endl; if (reproSys == 5 || reproSys == 6 || reproSys == 7) { cout << "Year PopSize nFemale " << endl; cout << "----------------------------" << endl; for (countY = 0; countY < (nYear+1); countY++) { cout << setw(8) << countY << " " << setw(8) << nPopSize[countY] << " " << setw(8) << nFemale[countY] << endl; } cout << endl; } else { cout << "Year PopSize " << endl; cout << "------------------" << endl; for (countY = 0; countY < (nYear+1); countY++) { cout << setw(8) << countY << " " << setw(8) << nPopSize[countY] << endl; } cout << endl; } } // end of displayVPopSizeInput function // displayGenotypeInput function void displayGenotypeInput(locusName *arrayInputLocus, individualId *arrayInputId, char ***arrayInput, int nSize, int nLoci, int allelePerLocus) { int countN; int countL; int countA; cout << "Imported data from input file: " << endl; cout << " Total number of individuals = " << nSize << endl; cout << " Total number of loci = " << nLoci << endl; cout << " Locus name list:"; for ( countL = 0; countL < nLoci; countL++) { cout << " " << arrayInputLocus[countL]; } cout << endl; cout << " Genotypic data:" << endl; for (countN = 0; countN < nSize; countN++) { cout << " " << arrayInputId[countN]; for ( countL = 0; countL < nLoci; countL++) { for ( countA = 0; countA < allelePerLocus; countA++) { cout << " " << arrayInput[countN][countL][countA]; } } cout << endl; } cout << endl; } // end of displayGenotypeInput function // checkAlleleFreqInput function void checkAlleleFreqInput(double *inputAlleleFreq, int nAllele) { int countA; double sumAlleleFreq = 0.0; // check for negative value for (countA = 0; countA < nAllele; countA++) { if (inputAlleleFreq[countA] < 0) { cout << endl; cout << "*** Input allele frequency error (negative value found!)."; exit(0); } } // check for allele frequency sum for (countA = 0; countA < nAllele; countA++) { sumAlleleFreq += inputAlleleFreq[countA]; } if (sumAlleleFreq <= 0.9999 || sumAlleleFreq >= 1.0001) { cout << endl; cout << "*** Input allele frequency error (not sum to 1!)."; exit(0); } } // end of checkAlleleFreqInput function // checkVPopSsizeInput function void checkVPopSsizeInput(int *nPopSize, int *nFemale, int reproSys, int nYear) { int countY; for (countY = 0; countY < (nYear+1); countY++) { switch (reproSys) { case 1: // asexual reproduction case 2: // Monoecious - complete selfing case 3: // Monoecious - random mating with selfing if (nPopSize[countY] < 1) { cout << "*** Population size setting error" << endl << "*** Population size in year" << countY << " < 1" << endl << "*** Population extinction occurred." << endl; exit(0); } break; case 4: // Monoecious - random mating without selfing if (nPopSize[countY] < 2) { cout << "*** Population size setting error" << endl << "*** Population size in year" << countY << " < 2" << endl << "*** Population extinction occurred." << endl; exit(0); } break; case 5: // Diecious - random mating case 6: // Diecious - single reproducing male case 7: // Diecious - single reproducing pair if (nPopSize[countY] < 2) // check nPopSize { cout << "*** Population size setting error" << endl << "*** Population size in year" << countY << " < 2" << endl << "*** Population extinction occurred." << endl; exit(0); } if (nFemale[countY] < 1) // check nFemale { cout << "*** Population size setting error" << endl << "*** nFemale in year" << countY << " < 1" << endl << "*** Population extinction occurred." << endl; exit(0); } if ((nPopSize[countY] - nFemale[countY]) < 1) // check nMale { cout << "*** Population size setting error" << endl << "*** nMale in year" << countY << " < 1" << endl << "*** Population extinction occurred." << endl; exit(0); } break; default: cout << endl; cout << "*** Reproductive system setting error" << endl; exit(0); } // end of reproSys switch } // end of countY loop } // end of checkVPopSsizeInput function // writeOutputHeader function void writeOutputHeader(ofstream &outputFile) { cout << "Writing header to output file..." << endl; cout << endl; // get time time_t rawtime; struct tm * timeinfo; time ( &rawtime ); timeinfo = localtime ( &rawtime ); // write header outputFile << "*** Start of output header ***" << endl; outputFile << "Program name: " << programName << endl; outputFile << "Version number: " << versionNumber << endl; outputFile << "Release date: " << versionDate << endl; outputFile << endl; outputFile << "Contact information" << endl; outputFile << " Author: " << author << endl; outputFile << " email: " << email << endl; outputFile << " Web: " << web << endl; outputFile << endl; outputFile << "Simulation date and time: " << asctime(timeinfo); outputFile << "*** End of output header ***" << endl; outputFile << endl; } // end of writeOutputHeader function // writeSLRawGenotypeHeader function void writeSLRawGenotypeHeader(ofstream &rawGenotypeOutputFile) { // get time time_t rawtime; struct tm * timeinfo; time ( &rawtime ); timeinfo = localtime ( &rawtime ); // write header rawGenotypeOutputFile << programName << " genotype raw data " << asctime(timeinfo); rawGenotypeOutputFile << "Loc1" << endl; } // end of writeSLRawGenotypeHeader function // writeMLRawGenotypeHeader function void writeMLRawGenotypeHeader(ofstream &rawGenotypeOutputFile, locusName *arrayInputLocus, int nLoci) { int countL; // get time time_t rawtime; struct tm * timeinfo; time ( &rawtime ); timeinfo = localtime ( &rawtime ); // write header rawGenotypeOutputFile << programName << " genotype raw data " << asctime(timeinfo); for (countL = 0; countL < nLoci; countL++) { rawGenotypeOutputFile << arrayInputLocus[countL] << endl; } } // end of writeMLRawGenotypeHeader function // writeVPopSizeInput function void writeVPopSizeInput(ofstream &outputFile, int *nPopSize, int *nFemale, int reproSys, int nYear) { int countY; outputFile << "*** Start of population size setting data from input file ***" << endl; outputFile << " Number of years = " << nYear << endl; outputFile << endl; if (reproSys == 5 || reproSys == 6 || reproSys == 7) { outputFile << "Year PopSize nFemale " << endl; outputFile << "----------------------------" << endl; for (countY = 0; countY < (nYear+1); countY++) { outputFile << setw(8) << countY << " " << setw(8) << nPopSize[countY] << " " << setw(8) << nFemale[countY] << endl; } outputFile << endl; } else { outputFile << "Year PopSize " << endl; outputFile << "------------------" << endl; for (countY = 0; countY < (nYear+1); countY++) { outputFile << setw(8) << countY << " " << setw(8) << nPopSize[countY] << endl; } outputFile << endl; } outputFile << "*** End of population size setting data from input file ***" << endl; } // end of writeVPopSizeInput function // writeAlleleFreq function // write allele freq data calculatted from input file to output file void writeAlleleFreq(ofstream &outputFile, int nLoci, locusName *arrayInputLocus, float **arrayAlleleFreq) { // declare variable // loop counting int countL; int countA; int countX; // output writing int boolPrint; // inputOA int *inputOA; inputOA = new int [nLoci]; if (inputOA == NULL) memoryError(); // inputEA float *inputEA; inputEA = new float [nLoci]; if (inputEA == NULL) memoryError(); // calculation float sumPi2 = 0; // sum (Pi)2 for calcaulating effective number of alleles // display cout << "Writing allele frequcy to output file..." << endl; // count number of alleles in each locus (nOA) for (countL = 0; countL < nLoci; countL++) { inputOA[countL] = 0; for (countA = 0; countA < maxAllele; countA++) { if (arrayAlleleFreq[countL][countA] > 0) inputOA[countL]++; } } // calculate nEA for (countL = 0; countL < nLoci; countL++) { sumPi2 = 0; for (countA = 0; countA < maxAllele; countA++) { sumPi2 += (arrayAlleleFreq[countL][countA] * arrayAlleleFreq[countL][countA]); } inputEA[countL] = (1 / sumPi2); } // write allele freq data to output file outputFile << "*** Start of allele frequency data from input file ***" << endl; outputFile << "------------"; for (countL = 0; countL < nLoci; countL++) { for (countX = 0; countX < locusNameLength; countX++) { outputFile << '-'; } } outputFile << endl; outputFile << "Allele/Locus"; for (countL = 0; countL < nLoci; countL++) { outputFile << setw(locusNameLength) << arrayInputLocus[countL]; } outputFile << endl; outputFile << "------------"; for (countL = 0; countL < nLoci; countL++) { for (countX = 0; countX < locusNameLength; countX++) { outputFile << '-'; } } outputFile << endl; for (countA = 0; countA < maxAllele; countA++) { // reset boolPrint boolPrint = 0; // check for allele presence in any locus for (countL = 0; countL < nLoci; countL++) { if (arrayAlleleFreq[countL][countA] != 0) { boolPrint = 1; } } if (boolPrint == 1) { outputFile << arrayAlleleId[countA] << " "; for (countL = 0; countL < nLoci; countL++) { if (arrayAlleleFreq[countL][countA] != 0) { outputFile << setw(locusNameLength) << setprecision(4) << setiosflags(ios::fixed | ios::showpoint) << arrayAlleleFreq[countL][countA]; } else { outputFile << setw(locusNameLength) << "0"; } } outputFile << endl; } } // write --- outputFile << "------------"; for (countL = 0; countL < nLoci; countL++) { for (countX = 0; countX < locusNameLength; countX++) { outputFile << '-'; } } outputFile << endl; // write nOA outputFile << "nObsAllele "; for (countL = 0; countL < nLoci; countL++) { outputFile << setw(locusNameLength) << inputOA[countL]; } outputFile << endl; // write --- outputFile << "------------"; for (countL = 0; countL < nLoci; countL++) { for (countX = 0; countX < locusNameLength; countX++) { outputFile << '-'; } } outputFile << endl; // write nEA outputFile << "nEffAllele "; for (countL = 0; countL < nLoci; countL++) { outputFile << setw(locusNameLength) << setprecision(4) << setiosflags(ios::fixed | ios::showpoint) << inputEA[countL]; } outputFile << endl; // write --- outputFile << "------------"; for (countL = 0; countL < nLoci; countL++) { for (countX = 0; countX < locusNameLength; countX++) { outputFile << '-'; } } outputFile << endl; outputFile << "*** End of allele frequency data from input file ***" << endl; outputFile << endl; } // end of writeAlleleFreq function // writeCPSSimSetting function void writeCPSSimSetting (ofstream &outputFile, int degreeOverlap, int reproSys, int longevity, int ageMaturation, int nPB, int nB, int nFemale, int nMale, int nYear, int nIteration) { // display message cout << "Writing simulation settings to output file..." << endl; cout << endl; outputFile << "*** Start of simulation settings ***" << endl; outputFile << "Simulation settings:" << endl; // call write functions writeDegreeOverlap(outputFile, degreeOverlap); writeReprosys(outputFile, reproSys); writeLongevity(outputFile, longevity); writeAgeMaturation(outputFile, ageMaturation); writeNPB(outputFile, nPB); writeNB(outputFile, nB); if ( reproSys == 5 || reproSys == 6 || reproSys == 7) { writeSexRatio(outputFile, nFemale, nMale); } writeNYear(outputFile, nYear); writeNIteration(outputFile, nIteration); outputFile << "*** End of simulation settings ***" << endl; } // end of writeCPSSimSetting function // writeVPSSimSetting function void writeVPSSimSetting (ofstream &outputFile, int degreeOverlap, int reproSys, int longevity, int ageMaturation, int nYear, int nIteration) { // display message cout << "Writing simulation settings to output file..." << endl; cout << endl; outputFile << "*** Start of simulation settings ***" << endl; outputFile << "Simulation settings:" << endl; // call write functions writeDegreeOverlap(outputFile, degreeOverlap); writeReprosys(outputFile, reproSys); writeLongevity(outputFile, longevity); writeAgeMaturation(outputFile, ageMaturation); writeNYear(outputFile, nYear); writeNIteration(outputFile, nIteration); outputFile << "*** End of simulation settings ***" << endl; } // end of writeVPSSimSetting function // writeDegreeOverlap function void writeDegreeOverlap(ofstream &outputFile, int degreeOverlap) { outputFile << " Generation overlap setting = " << degreeOverlap << endl; outputFile << " (0) = Minimum overlap (discrete generations)" << endl; outputFile << " (100) = Maximum overlap (random starting age for all individuals)" << endl; outputFile << endl; } // end of writeDegreeOverlap function // writeReprosys function void writeReprosys(ofstream &outputFile, int reproSys) { outputFile << " Reproductive system setting = " << reproSys << endl; outputFile << " (1) = Asexual reproduction" << endl; outputFile << " (2) = Monoecy with complete selfing" << endl; outputFile << " (3) = Monoecy with random mating (with selfing)" << endl; outputFile << " (4) = Monoecy with random mating (without selfing)" << endl; outputFile << " (5) = Dioecy with random mating" << endl; outputFile << " (6) = Dioecy with single reproducing male each year" << endl; outputFile << " (7) = Dioecy with single reproducing pair each year" << endl; outputFile << endl; } // end of writeReprosys function // writeLongevity function void writeLongevity(ofstream &outputFile, int longevity) { outputFile << " Expected longevity of the organism = " << longevity << " year(s)" < 0) nOA++; } // check if fixed if (nOA == 1) { nEA = 1; nHo = 0; nHe = 0; nF = 1; } else { // allele freq for (countA = 0; countA < nAllele; countA++) { arrayAlleleFreq[countA] = float (arrayAlleleCount[countA]) / (nPopSize * allelePerLocus); } // effective number of alleles sumPi2 = 0; // reset for (countA = 0; countA < nAllele; countA++) { sumPi2 += (arrayAlleleFreq[countA] * arrayAlleleFreq[countA]); } nEA = 1 / sumPi2; // Ho nHeteroFound = 0; // reset for (countN = 0; countN < nPopSize; countN++) { if (arrayGenotype[countN][0] != arrayGenotype[countN][1]) nHeteroFound++; } nHo = float (nHeteroFound) / nPopSize; // He nHe = 1; // reset for (countA = 0; countA < nAllele; countA++) { nHe -= (arrayAlleleFreq[countA] * arrayAlleleFreq[countA]); } if (nHe < 0) { nHe = -nHe; } // F nF = (nHe - nHo) / nHe; } } // end of statDiploidSingleLocus function // statDiploidMultiLocus function // receive genotype array from simulation module // receive nPopSize (# of individuals) and nLoci (# of Loci) // return stat (nOA, nEA, nHo, nHe, nF) void statDiploidMultiLocus(int ***arrayGenotype, int allelePerLocus, int nPopSize, int nLoci, int *arrayOA, float *arrayEA, float *arrayHo, float *arrayHe, float *arrayF) { // declare variables // loop counting int countN; int countL; int countA; // calculation int nHeteroFound; // observed number of heterozygous individuals float sumPi2 = 0; // sum (Pi)2 for calcaulating effective number of alleles // initiate array // allele count array [maxAllele] int *arrayAlleleCount; arrayAlleleCount = new int [maxAllele]; if (arrayAlleleCount == NULL) memoryError(); // allele freq array [maxAllele] float *arrayAlleleFreq; arrayAlleleFreq = new float [maxAllele]; if (arrayAlleleFreq == NULL) memoryError(); // stat calculation for (countL = 0; countL < nLoci; countL++) { // allele count for ( countA = 0; countA < maxAllele; countA++) { arrayAlleleCount[countA] = 0; } for (countN = 0; countN < nPopSize; countN++) { for (countA = 0; countA < allelePerLocus; countA++) { arrayAlleleCount[(arrayGenotype[countN][countL][countA])]++; } } // OA arrayOA[countL] = 0; // reset for (countA = 0; countA < maxAllele; countA++) { if (arrayAlleleCount[countA] > 0) arrayOA[countL]++; } // check if fixed if (arrayOA[countL] == 1) { arrayEA[countL] = 1; arrayHo[countL] = 0; arrayHe[countL] = 0; arrayF[countL] = 1; } else { // allele freq for (countA = 0; countA < maxAllele; countA++) { arrayAlleleFreq[countA] = float (arrayAlleleCount[countA]) / (nPopSize * allelePerLocus); } // effective number of alleles sumPi2 = 0; // reset for (countA = 0; countA < maxAllele; countA++) { sumPi2 += (arrayAlleleFreq[countA] * arrayAlleleFreq[countA]); } arrayEA[countL] = 1 / sumPi2; // Ho nHeteroFound = 0; // reset for (countN = 0; countN < nPopSize; countN++) { if (arrayGenotype[countN][countL][0] != arrayGenotype[countN][countL][1]) nHeteroFound++; } arrayHo[countL] = float (nHeteroFound) / nPopSize; // He arrayHe[countL] = 1; // reset for (countA = 0; countA < maxAllele; countA++) { arrayHe[countL] -= (arrayAlleleFreq[countA] * arrayAlleleFreq[countA]); } if (arrayHe[countL] < 0) { arrayHe[countL] = -arrayHe[countL]; } // F arrayF[countL] = (arrayHe[countL] - arrayHo[countL]) / arrayHe[countL]; } } // end of countL loop } // end of statDiploidMultiLocus function // outputDiploidSingleLocus function (int) // receive raw data array (int), nYear, nIteration from simulation module // calculate and return mean, StdD, StdE, and % array void outputDiploidSingleLocus(int **arrayRawdata, int nYear, int nIteration, float *arrayMean, float *arrayStdD, float *arrayStdE, float *arrayPercent) { // declare variables // loop counting int countY; int countI; // calculation int total; float sumS; // calculation // calculate mean for (countY = 0; countY < (nYear + 1); countY++) { total = 0; for (countI = 0; countI < nIteration; countI++) { total += arrayRawdata[countY][countI]; } arrayMean[countY] = ( float (total) / nIteration); } // calculate stdD, stdE, % for (countY = 0; countY < (nYear + 1); countY++) { sumS = 0; for (countI = 0; countI < nIteration; countI++) { sumS += (float (arrayRawdata[countY][countI]) - arrayMean[countY]) * (float (arrayRawdata[countY][countI]) - arrayMean[countY]); } arrayStdD[countY] = (sqrt((sumS/(nIteration - 1)))); arrayStdE[countY] = (arrayStdD[countY] / sqrt((float(nIteration)))); arrayPercent[countY] = ((arrayMean[countY] / arrayMean[0]) * 100); } } // end of outputDiploidSingleLocus function (int) // outputDiploidSingleLocus function (float) // receive raw data array (float), nYear, nIteration from simulation module // calculate and return mean, StdD, StdE, and % array void outputDiploidSingleLocus(float **arrayRawdata, int nYear, int nIteration, float *arrayMean, float *arrayStdD, float *arrayStdE, float *arrayPercent) { // declare variables // loop counting int countY; int countI; // calculation float total; float sumS; // calculation // calculate mean for (countY = 0; countY < (nYear + 1); countY++) { total = 0; for (countI = 0; countI < nIteration; countI++) { total += arrayRawdata[countY][countI]; } arrayMean[countY] = (total / nIteration); } // calculate stdD, stdE, % for (countY = 0; countY < (nYear + 1); countY++) { sumS = 0; for (countI = 0; countI < nIteration; countI++) { sumS += (arrayRawdata[countY][countI] - arrayMean[countY]) * (arrayRawdata[countY][countI] - arrayMean[countY]); } arrayStdD[countY] = (sqrt((sumS/(nIteration - 1)))); arrayStdE[countY] = (arrayStdD[countY] / sqrt(float(nIteration))); arrayPercent[countY] = ((arrayMean[countY] / arrayMean[0]) * 100); } } // end of outputDiploidSingleLocus function (float) // outputDiploidSingleLocus function (int) // receive raw data array (int), nLoci, nYear, nIteration from simulation module // calculate and return mean, Av, StdD, StdE, and % array void outputDiploidMultiLocus(int ***arrayRawdata, int nLoci, int nYear, int nIteration, float **arrayMean, float *arrayAv, float *arrayStdD, float *arrayStdE, float *arrayPercent) { // declare variables // loop counting int countL; int countY; int countI; // calculation int intTotal; float floatTotal; float sumS; // calculation // calculate mean for (countL = 0; countL < nLoci; countL++) { for (countY = 0; countY < (nYear + 1); countY++) { intTotal = 0; for (countI = 0; countI < nIteration; countI++) { intTotal += arrayRawdata[countL][countY][countI]; } arrayMean[countL][countY] = ( float (intTotal) / nIteration); } } // calculate Av for (countY = 0; countY < (nYear + 1); countY++) { floatTotal = 0; for (countL = 0; countL < nLoci; countL++) { floatTotal += arrayMean[countL][countY]; } arrayAv[countY] = (floatTotal / nLoci); } // calculate stdD, stdE, % for (countY = 0; countY < (nYear + 1); countY++) { sumS = 0; for (countL = 0; countL < nLoci; countL++) { sumS += ((arrayMean[countL][countY] - arrayAv[countY]) * (arrayMean[countL][countY] - arrayAv[countY])); } arrayStdD[countY] = (sqrt((sumS / (nLoci - 1)))); arrayStdE[countY] = (arrayStdD[countY] / sqrt(float(nLoci))); arrayPercent[countY] = ((arrayAv[countY] / arrayAv[0]) * 100); } } // end of outputDiploidSingleLocus function (int) // outputDiploidSingleLocus function (float) // receive raw data array (float), nLoci, nYear, nIteration from simulation module // calculate and return mean, Av, StdD, StdE, and % array void outputDiploidMultiLocus(float ***arrayRawdata, int nLoci, int nYear, int nIteration, float **arrayMean, float *arrayAv, float *arrayStdD, float *arrayStdE, float *arrayPercent) { // declare variables // loop counting int countL; int countY; int countI; // calculation float floatTotal; float sumS; // calculation // calculate mean for (countL = 0; countL < nLoci; countL++) { for (countY = 0; countY < (nYear + 1); countY++) { floatTotal = 0; for (countI = 0; countI < nIteration; countI++) { floatTotal += arrayRawdata[countL][countY][countI]; } arrayMean[countL][countY] = (floatTotal / nIteration); } } // calculate Av for (countY = 0; countY < (nYear + 1); countY++) { floatTotal = 0; for (countL = 0; countL < nLoci; countL++) { floatTotal += arrayMean[countL][countY]; } arrayAv[countY] = (floatTotal / nLoci); } // calculate stdD, stdE, % for (countY = 0; countY < (nYear + 1); countY++) { sumS = 0; for (countL = 0; countL < nLoci; countL++) { sumS += ((arrayMean[countL][countY] - arrayAv[countY]) * (arrayMean[countL][countY] - arrayAv[countY])); } arrayStdD[countY] = (sqrt((sumS / (nLoci - 1)))); arrayStdE[countY] = (arrayStdD[countY] / sqrt(float(nLoci))); arrayPercent[countY] = ((arrayAv[countY] / arrayAv[0]) * 100); } } // end of outputDiploidSingleLocus function (float) // openInput function // open input file and get input file name void openInput(ifstream &inputFile, fileName &inputFileName) { // prompt for input filename cout << endl; cout << "Input file must be in the same directory as this program," << endl; cout << "The limitation of file name length is " << (fileNameLength - 1) << " characters." << endl; cout << "Please enter input filename: "; cin >> inputFileName; // open input file inputFile.open(inputFileName); // check input file if ( inputFile == 0 ) { cout << endl; cout << "*** Error opening input file." << endl; exit(0); } else cout << endl; cout << "Input file (" << inputFileName << ") opened successfully." << endl; } // end of openInput // openPopSizeInput function // open input file and get input file name void openPopSizeInput(ifstream &popSizeInput, fileName &popSizeInputName) { // prompt for input filename cout << endl; cout << "Population size setting file must be in the same directory as this program," << endl; cout << "The limitation of file name length is " << (fileNameLength - 1) << " characters." << endl; cout << "Please enter population size setting filename: "; cin >> popSizeInputName; // open input file popSizeInput.open(popSizeInputName); // check input file if ( popSizeInput == 0 ) { cout << endl; cout << "*** Error opening population size setting file." << endl; exit(0); } else cout << endl; cout << "Population size setting file (" << popSizeInputName << ") opened successfully." << endl; } // end of openPopSizeInput // openOutput function // open output file and get output file name void openOutput(ofstream &outputFile, fileName &outputFileName) { // prompt for output filename cout << endl; cout << "The limitation of file name length is " << (fileNameLength - 1) << " characters." << endl; cout << "Please enter output filename: "; cin >> outputFileName; // open outputFile outputFile.open(outputFileName); if ( outputFile == 0 ) { cout << endl; cout << "*** Error opening output file." << endl; exit(0); } else cout << endl; cout << "Output file (" << outputFileName << ") opened successfully." << endl; } // end of openOutput // openRawGenotypeOutput function void openRawGenotypeOutput(ofstream &rawGenotypeOutputFile, fileName &rawGenotypeOutputName) { // prompt for output filename cout << endl; cout << "The limitation of file name length is " << (fileNameLength - 1) << " characters." << endl; cout << "Please enter raw genotype data output filename: "; cin >> rawGenotypeOutputName; // open outputFile rawGenotypeOutputFile.open(rawGenotypeOutputName); if ( rawGenotypeOutputFile == 0 ) { cout << endl; cout << "*** Error opening output file." << endl; exit(0); } else cout << endl; cout << "Output file (" << rawGenotypeOutputName << ") opened successfully." << endl; } // end of openRawGenotypeOutput function // closeInputFile function // close input file void closeInputFile(ifstream &inputFile) { cout << endl; inputFile.close(); } // closeOutputFile function // close output file void closeOutputFile(ofstream &outputFile) { cout << endl; outputFile.close(); } // simDiploidSingleLocusCPS function // perform diploid singlelocus constant population size simulation void simDiploidSingleLocusCPS(void) { // constants const int allelePerLocus = 2; // number of alleles per locus // declare loop-counting variables int countN = 0; // count individual int countNTemp = 0; // count individual int countA = 0; // count allele int countAPool = 0; // count allele pool int countY = 0; // count year int countI = 0; // count iteration // file operation // input file - allele freq setting // declare fileName fileName inputFileName; // input - allele freq // declare input file ifstream alleleFreqInput; // call openInput function openInput(alleleFreqInput, inputFileName); // output file - simulation summary // declare fileName fileName outputFileName; // output - simulation summary // declare output file ofstream outputFile; // call openOutput function openOutput(outputFile, outputFileName); // call writeOutputHeader function // write header to output file writeOutputHeader(outputFile); // output file - raw genotype data // declare fileName fileName rawGenotypeOutputName; // declare raw data output file ofstream rawGenotypeOutputFile; // read input file // declare variables int nAllele = 0; // observed number of alleles in the pre bottleneck population // display cout << endl; cout << "Reading input file..." << endl; // read nAllele nAllele = readInputFileNAllele(alleleFreqInput); // initiate input array (input allele frequency) double *inputAlleleFreq; // array for input allele frequency [nAllele] inputAlleleFreq = new double[nAllele]; if (inputAlleleFreq == NULL) memoryError(); // read input allele freq readInputFileAlleleFreq(alleleFreqInput, inputAlleleFreq, nAllele); // display imported input data displayAlleleFreqInput(inputAlleleFreq, nAllele); // check input allele freq data checkAlleleFreqInput(inputAlleleFreq, nAllele); // get simulation settings // declare variables int degreeOverlap; // generation overlapping setting int reproSys; // reproductive system option int longevity; // longevity of the organism int ageMaturation; // reproductive maturation age of the organism int nPB; // pre-bottleneck population size int nB; // bottlneck population size // number of female/male for diecious models int nFemale; int nMale; int nYear; // duration of bottleneck int nIteration; // number of iterations of the run int boolOutputRawGenotype; // bool variable for output raw genotype data // call getCPSSimSetting function getCPSSimSetting(degreeOverlap, reproSys, longevity, ageMaturation, nPB, nB, nFemale, nMale, nYear, nIteration, boolOutputRawGenotype); // check boolOutputRawGenotype if (boolOutputRawGenotype == 1) { // open raw data output file // call openRawGenotypeOutput function openRawGenotypeOutput(rawGenotypeOutputFile, rawGenotypeOutputName); // write header writeSLRawGenotypeHeader(rawGenotypeOutputFile); } // declare simulation variables int parent1 = 0; // id of parent1 int parent2 = 0; // id of parent2 int nSurvivedFemale = 0; // number of survived females int nSurvivedMale = 0; // number of survived males int nSurvivedTotal = 0; // number of survivors int nAdultFemale = 0; // number of adult females int nAdultMale = 0; // number of adult Males int nAdultTotal = 0; // number of adults // initiate data arrays // pre-bottleneck population genotype [nPB][allelePerLocus] int **arrayPBP; arrayPBP = new int* [nPB]; if (arrayPBP == NULL) memoryError(); for (countN = 0; countN < nPB; countN++) { arrayPBP[countN] = new int[allelePerLocus]; if (arrayPBP == NULL) memoryError(); } // bottleneck population genotype [nB][allelePerLocus] int **arrayBP; arrayBP = new int* [nB]; if (arrayBP == NULL) memoryError(); for (countN = 0; countN < nB; countN++) { arrayBP[countN] = new int[allelePerLocus]; if (arrayBP == NULL) memoryError(); } // breeding population genotype [nB][allelePerLocus] int **arrayParent; arrayParent = new int* [nB]; if (arrayParent == NULL) memoryError(); for (countN = 0; countN < nB; countN++) { arrayParent[countN] = new int[allelePerLocus]; if (arrayParent == NULL) memoryError(); } // new genotype array int *arrayNewGenotype; arrayNewGenotype = new int [allelePerLocus]; if (arrayNewGenotype == NULL) memoryError(); // age [nB] int *arrayAge; arrayAge = new int [nB]; if (arrayAge == NULL) memoryError(); // allele count array [nAllele] int *arrayAlleleCount; arrayAlleleCount = new int [nAllele]; if (arrayAlleleCount == NULL) memoryError(); // declare variable int nOA = 0; // observed number of alleles float nEA = 0; // effective number of alleles float nHo = 0; // observed heterozygosity float nHe = 0; // expected heterozygosity float nF = 0; // fixation index // initiate raw data array // OA [(nYear+1)][nIteration] int **arrayRawOA; arrayRawOA = new int* [(nYear+1)]; if (arrayRawOA == NULL) memoryError(); for (countY = 0; countY < (nYear + 1); countY++) { arrayRawOA[countY] = new int[nIteration]; if (arrayRawOA == NULL) memoryError(); } // EA [(nYear+1)][nIteration] float **arrayRawEA; arrayRawEA = new float* [(nYear+1)]; if (arrayRawEA == NULL) memoryError(); for (countY = 0; countY < (nYear + 1); countY++) { arrayRawEA[countY] = new float[nIteration]; if (arrayRawEA == NULL) memoryError(); } // Ho [(nYear+1)][nIteration] float **arrayRawHo; arrayRawHo = new float* [(nYear+1)]; if (arrayRawHo == NULL) memoryError(); for (countY = 0; countY < (nYear + 1); countY++) { arrayRawHo[countY] = new float[nIteration]; if (arrayRawHo == NULL) memoryError(); } // He [(nYear+1)][nIteration] float **arrayRawHe; arrayRawHe = new float* [(nYear+1)]; if (arrayRawHe == NULL) memoryError(); for (countY = 0; countY < (nYear + 1); countY++) { arrayRawHe[countY] = new float[nIteration]; if (arrayRawHe == NULL) memoryError(); } // F [(nYear+1)][nIteration] float **arrayRawF; arrayRawF = new float* [(nYear+1)]; if (arrayRawF == NULL) memoryError(); for (countY = 0; countY < (nYear + 1); countY++) { arrayRawF[countY] = new float[nIteration]; if (arrayRawF == NULL) memoryError(); } // vector for list // allele pool [(nPB * allelePerLocus)] // allele pool of the pre-bottleneck population vector vectAllelePool; // get allele count getDiSLAlleleCount(arrayAlleleCount, inputAlleleFreq, nAllele, nPB, allelePerLocus); // fill allele pool vector for (countA = 0; countA < nAllele; countA++) { for (countAPool = 0; countAPool < (arrayAlleleCount[countA]); countAPool++) { vectAllelePool.push_back(countA); } } // pre bottleneck population list [nPB] // list of all individuals in the pre bottleneck population // for selecting individuals survived bottleneck vector vectPBPList; for (countN = 0; countN < nPB; countN++) vectPBPList.push_back(countN); // start of iteration loop for (countI = 0; countI < nIteration; countI++) { // generation pre-bottleneck population // random shuffle allele pool random_shuffle(vectAllelePool.begin(), vectAllelePool.end()); // fill genotype array countAPool = 0; // reset allele pool counting variable for (countN = 0; countN < nPB; countN++) { for (countA = 0; countA < allelePerLocus; countA++) { arrayPBP[countN][countA] = vectAllelePool[countAPool]; countAPool++; } } // stat of pre bottleneck population // call statDiploidSingleLocus function statDiploidSingleLocus(arrayPBP, allelePerLocus, nPB, nAllele, nOA, nEA, nHo, nHe, nF); // write to raw data array arrayRawOA[0][countI] = nOA; arrayRawEA[0][countI] = nEA; arrayRawHo[0][countI] = nHo; arrayRawHe[0][countI] = nHe; arrayRawF[0][countI] = nF; // generate bottleneck population random_shuffle(vectPBPList.begin(), vectPBPList.end()); // survivor list for (countN = 0; countN < nB; countN++) { // genotype for (countA = 0; countA < allelePerLocus; countA++) { arrayBP[countN][countA] = arrayPBP[(vectPBPList[countN])][countA]; } // age // when degreeOverlap = 0, if-test will always return true -> all individual start at age 0 // when degreeOverlap = 100, if-test will always return fales -> all individual start with random age value if (degreeOverlap < getRandomNumber(1, 100)) { arrayAge[countN] = 0; } else { arrayAge[countN] = getRandomNumber(0, (longevity - 1)); } } // start of year loop // each year loop has 2 major steps // 1. from previous year end to current year begin // generate survived population from the previous year // 2. from current year begin to current year end // generate breeding population of the current year // check age/replace old individuals reached longevity limit // geneerate new individuals if population growth occurred // // individual death events (old genotypes eliminated or replaced by new ones) can happen at both steps // death in step 1 is due to population decline (random death) // death in step 2 is due to reaching longevity limit (check age) // individual birth events always happen at step 2 // possible birth event1: old individuals reached longevity limit // possible birth event2: population growth occurred for (countY = 0; countY < nYear; countY++) { // check if asexual/monoecious or dioecious // if asexual/monoecious --> treat the population as a whole // else (dioecious) --> treat female and male separately if (reproSys >= 1 && reproSys <= 4) // asexual/monoecious, treat the population aas a whole { // step1: from previous year end to curent year begin // in constant population size module all individuals survived // do nothing to genotype array nSurvivedTotal = nB; // age increased by 1 for (countN = 0; countN < nSurvivedTotal; countN++) { arrayAge[countN]++; } // step2: from current year begin to current year end // generate breeding population // check individual age // if age reached longevity limit -> generate new individuals // generate breeding population // check if reached reproductive maturation, copy genotype to arrayParent nAdultTotal = 0; // reset number of adults for (countN = 0; countN < nSurvivedTotal; countN++) { if (arrayAge[countN] >= ageMaturation) { arrayParent[nAdultTotal][0] = arrayBP[countN][0]; arrayParent[nAdultTotal][1] = arrayBP[countN][1]; nAdultTotal++; } } // check survivor age // reset age and replace with new genotype if reached limit // else do nothing to genotype and age array for (countN = 0; countN < nSurvivedTotal; countN++) { // check if reached longevity limit, generate new individual if (arrayAge[countN] >= longevity) { // reset age arrayAge[countN] = 0; // generate new genotype // total nAdultTotal breeders getNewDiSLGenotype(arrayParent, arrayNewGenotype, reproSys, nAdultTotal, nAdultFemale, parent1, parent2); // copy genotype arrayBP[countN][0] = arrayNewGenotype[0]; arrayBP[countN][1] = arrayNewGenotype[1]; } // end of check age if // else do nothing to genotype and age array } // end of check age for loop } // end of asexual/monoecious if else // dioecious reproSys, treat female/male separately { // step1: from previous year end to curent year begin // in constant population size module all individuals survived // do nothing to genotype array nSurvivedFemale = nFemale; nSurvivedMale = nMale; nSurvivedTotal = nB; // age increased by 1 for (countN = 0; countN < nSurvivedTotal; countN++) { arrayAge[countN]++; } // step2: from current year begin to current year end // generate breeding population // check individual age // if age reached longevity limit -> generate new individuals // generate breeding population // check if reached reproductive maturation, copy genotype to arrayParent // generate breeding population nAdultFemale = 0; for (countN = 0; countN < nSurvivedFemale; countN++) { if (arrayAge[countN] >= ageMaturation) { arrayParent[nAdultFemale][0] = arrayBP[countN][0]; arrayParent[nAdultFemale][1] = arrayBP[countN][1]; nAdultFemale++; } } nAdultMale = 0; for (countN = 0; countN < nSurvivedMale; countN++) { if (arrayAge[(countN + nSurvivedFemale)] >= ageMaturation) { arrayParent[(nAdultMale + nAdultFemale)][0] = arrayBP[(countN + nSurvivedFemale)][0]; arrayParent[(nAdultMale + nAdultFemale)][1] = arrayBP[(countN + nSurvivedFemale)][1]; nAdultMale++; } } nAdultTotal = (nAdultFemale + nAdultMale); // get parent id for diecious model with single reproducing male/pair // single reproducing male // get male id if (reproSys == 6) { parent2 = getRandomNumber(nAdultFemale, (nAdultTotal - 1)); } // single reproducing pair // get both parent ids if (reproSys == 7) { parent1 = getRandomNumber(0, (nAdultFemale - 1)); parent2 = getRandomNumber(nAdultFemale, (nAdultTotal - 1)); } // check survivor age // reset age and replace with new genotype if reached limit // else do nothing to genotype and age array for (countN = 0; countN < nSurvivedTotal; countN++) { // check if reached longevity limit, generate new individual if (arrayAge[countN] >= longevity) { // reset age arrayAge[countN] = 0; // generate new genotype // total nAdultTotal breeders getNewDiSLGenotype(arrayParent, arrayNewGenotype, reproSys, nAdultTotal, nAdultFemale, parent1, parent2); // copy genotype arrayBP[countN][0] = arrayNewGenotype[0]; arrayBP[countN][1] = arrayNewGenotype[1]; } // end of check age if // else do nothing to genotype and age array } // end of check age for loop } // end of dioecious repro sys else // stat of bottleneck population // call statDiploidSingleLocus function statDiploidSingleLocus(arrayBP, allelePerLocus, nB, nAllele, nOA, nEA, nHo, nHe, nF); // write to data array arrayRawOA[(countY + 1)][countI] = nOA; arrayRawEA[(countY + 1)][countI] = nEA; arrayRawHo[(countY + 1)][countI] = nHo; arrayRawHe[(countY + 1)][countI] = nHe; arrayRawF[(countY + 1)][countI] = nF; } // end of year loop // check boolOutputRawGenotype if (boolOutputRawGenotype == 1) { // call writeDiSLRawGenotype function writeDiSLRawGenotype(rawGenotypeOutputFile, arrayBP, allelePerLocus, nB); } cout << '.'; // use dot to show iteration progress } // end of iteration loop // display message cout << endl; cout << "Simulation completed." << endl; cout << endl; cout << "Preparing output file..." << endl; cout << endl; // write simulation settings to output file // call writeCPSSimSetting function writeCPSSimSetting(outputFile, degreeOverlap, reproSys, longevity, ageMaturation, nPB, nB, nFemale, nMale, nYear, nIteration); // write summary output file writeDiSLSummary(outputFile, nYear, nIteration, arrayRawOA, arrayRawEA, arrayRawHo, arrayRawHe, arrayRawF); // call closeFile function // close input file closeInputFile(alleleFreqInput); // close output file closeOutputFile(outputFile); // check boolOutputRawGenotype if (boolOutputRawGenotype == 1) { closeOutputFile(rawGenotypeOutputFile); } } // end of simDiploidSingleLocusCPS function // simDiploidSingleLocusVPS function // perform diploid singlelocus variable population size simulation void simDiploidSingleLocusVPS(void) { // constants const int allelePerLocus = 2; // number of alleles per locus // declare loop-counting variables int countN = 0; // count individual int countNTemp = 0; // count individual int countA = 0; // count allele int countAPool = 0; // count allele pool int countY = 0; // count year int countI = 0; // count iteration // file operation // input file - allele freq setting // declare fileName fileName alleleFreqInputName; // input - allele freq // declare input file ifstream alleleFreqInput; // call openInput function openInput(alleleFreqInput, alleleFreqInputName); // input file - population size setting // declare fileName fileName popSizeInputName; // input - population size // declare input file ifstream popSizeInput; // output file - simulation summary // declare fileName fileName outputFileName; // output - simulation summary // declare output file ofstream outputFile; // call openOutput function openOutput(outputFile, outputFileName); // call writeOutputHeader function // write header to output file writeOutputHeader(outputFile); // output file - raw genotype data // declare fileName fileName rawGenotypeOutputName; // declare raw data output file ofstream rawGenotypeOutputFile; // read input file // declare variables int nAllele = 0; // observed number of alleles in the pre bottleneck population // display cout << endl; cout << "Reading allele frequency input file..." << endl; // read nAllele nAllele = readInputFileNAllele(alleleFreqInput); // initiate input array (input allele frequency) double *inputAlleleFreq; // array for input allele frequency [nAllele] inputAlleleFreq = new double[nAllele]; if (inputAlleleFreq == NULL) memoryError(); // read input allele freq readInputFileAlleleFreq(alleleFreqInput, inputAlleleFreq, nAllele); // display imported input data displayAlleleFreqInput(inputAlleleFreq, nAllele); // check input allele freq data checkAlleleFreqInput(inputAlleleFreq, nAllele); // get simulation settings // declare variables int degreeOverlap; // generation overlapping setting int reproSys; // reproductive system option int longevity; // longevity of the organism int ageMaturation; // reproductive maturation age of the organism int nIteration; // number of iterations of the run int boolOutputRawGenotype; // bool variable for output raw genotype data // call getVPSSimSetting function getVPSSimSetting(degreeOverlap, reproSys, longevity, ageMaturation, nIteration, boolOutputRawGenotype); // check boolOutputRawGenotype if (boolOutputRawGenotype == 1) { // open raw data output file // call openRawGenotypeOutput function openRawGenotypeOutput(rawGenotypeOutputFile, rawGenotypeOutputName); // write header writeSLRawGenotypeHeader(rawGenotypeOutputFile); } // get population size settings // call openPopSizeInput function openPopSizeInput(popSizeInput, popSizeInputName); // read population size setting input file // declare variables int nYear; // duration of bottleneck // display cout << endl; cout << "Reading population size setting input file..." << endl; // read nYear nYear = readInputFileNYear(popSizeInput); // initiate pop size input data array int *nPopSize; // array for popsize [(nYear+1)] nPopSize = new int [(nYear+1)]; if (nPopSize == NULL) memoryError(); int *nFemale; // array for nFemale [(nYear+1)] nFemale = new int [(nYear+1)]; if (nFemale == NULL) memoryError(); int *nMale; // array for nMale [(nYear+1)] nMale = new int [(nYear+1)]; if (nMale == NULL) memoryError(); // read population size input readInputFilePopSize(popSizeInput, nPopSize, nFemale, nMale, nYear, reproSys); // write imported population size input data writeVPopSizeInput(outputFile, nPopSize, nFemale, reproSys, nYear); // display imported population size input data displayVPopSizeInput(nPopSize, nFemale, reproSys, nYear); // check population size input checkVPopSsizeInput(nPopSize, nFemale, reproSys, nYear); // get max pop size int nMaxPopSize; // maximum population size in pop size input nMaxPopSize = getNMaxPopSize(nPopSize, nYear); // declare simulation variables int parent1 = 0; // id of parent1 int parent2 = 0; // id of parent2 int nSurvivedFemale = 0; // number of survived females int nSurvivedMale = 0; // number of survived males int nSurvivedTotal = 0; // number of survivors int nAdultFemale = 0; // number of adult females int nAdultMale = 0; // number of adult Males int nAdultTotal = 0; // number of adults // initiate data arrays // population genotype [nMaxPopSize][allelePerLocus] int **arrayGenotype; arrayGenotype = new int* [nMaxPopSize]; if (arrayGenotype == NULL) memoryError(); for (countN = 0; countN < nMaxPopSize; countN++) { arrayGenotype[countN] = new int[allelePerLocus]; if (arrayGenotype == NULL) memoryError(); } // population genotype (temp) [nMaxPopSize][allelePerLocus] int **arrayGenotypeTemp; arrayGenotypeTemp = new int* [nMaxPopSize]; if (arrayGenotypeTemp == NULL) memoryError(); for (countN = 0; countN < nMaxPopSize; countN++) { arrayGenotypeTemp[countN] = new int[allelePerLocus]; if (arrayGenotypeTemp == NULL) memoryError(); } // breeding population genotype [nMaxPopSize][allelePerLocus] int **arrayParent; arrayParent = new int* [nMaxPopSize]; if (arrayParent == NULL) memoryError(); for (countN = 0; countN < nMaxPopSize; countN++) { arrayParent[countN] = new int[allelePerLocus]; if (arrayParent == NULL) memoryError(); } // new genotype int *arrayNewGenotype; arrayNewGenotype = new int [allelePerLocus]; if (arrayNewGenotype == NULL) memoryError(); // age [nMaxPopSize] int *arrayAge; arrayAge = new int [nMaxPopSize]; if (arrayAge == NULL) memoryError(); // age (temp) [nMaxPopSize] int *arrayAgeTemp; arrayAgeTemp = new int [nMaxPopSize]; if (arrayAgeTemp == NULL) memoryError(); // allele count array [nAllele] int *arrayAlleleCount; arrayAlleleCount = new int [nAllele]; if (arrayAlleleCount == NULL) memoryError(); // declare variable int nOA = 0; // observed number of alleles float nEA = 0; // effective number of alleles float nHo = 0; // observed heterozygosity float nHe = 0; // expected heterozygosity float nF = 0; // fixation index // initiate raw data array // OA [(nYear+1)][nIteration] int **arrayRawOA; arrayRawOA = new int* [(nYear+1)]; if (arrayRawOA == NULL) memoryError(); for (countY = 0; countY < (nYear + 1); countY++) { arrayRawOA[countY] = new int[nIteration]; if (arrayRawOA == NULL) memoryError(); } // EA [(nYear+1)][nIteration] float **arrayRawEA; arrayRawEA = new float* [(nYear+1)]; if (arrayRawEA == NULL) memoryError(); for (countY = 0; countY < (nYear + 1); countY++) { arrayRawEA[countY] = new float[nIteration]; if (arrayRawEA == NULL) memoryError(); } // Ho [(nYear+1)][nIteration] float **arrayRawHo; arrayRawHo = new float* [(nYear+1)]; if (arrayRawHo == NULL) memoryError(); for (countY = 0; countY < (nYear + 1); countY++) { arrayRawHo[countY] = new float[nIteration]; if (arrayRawHo == NULL) memoryError(); } // He [(nYear+1)][nIteration] float **arrayRawHe; arrayRawHe = new float* [(nYear+1)]; if (arrayRawHe == NULL) memoryError(); for (countY = 0; countY < (nYear + 1); countY++) { arrayRawHe[countY] = new float[nIteration]; if (arrayRawHe == NULL) memoryError(); } // F [(nYear+1)][nIteration] float **arrayRawF; arrayRawF = new float* [(nYear+1)]; if (arrayRawF == NULL) memoryError(); for (countY = 0; countY < (nYear + 1); countY++) { arrayRawF[countY] = new float[nIteration]; if (arrayRawF == NULL) memoryError(); } // vector for list // allele pool [(nPopSize[0] * allelePerLocus)] // allele pool of the pre-bottleneck population vector vectAllelePool; // get allele count getDiSLAlleleCount(arrayAlleleCount, inputAlleleFreq, nAllele, nPopSize[0], allelePerLocus); // fill allele pool vector for (countA = 0; countA < nAllele; countA++) { for (countAPool = 0; countAPool < (arrayAlleleCount[countA]); countAPool++) { vectAllelePool.push_back(countA); } } // population list // list of all individuals in the population // for selecting surviving individuals vector vectPopList; // start of iteration loop for (countI = 0; countI < nIteration; countI++) { // generation year0 population // random shuffle allele pool random_shuffle(vectAllelePool.begin(), vectAllelePool.end()); // fill genotype array (previous year end genotype array) countAPool = 0; // reset allele pool counting variable for (countN = 0; countN < nPopSize[0]; countN++) { for (countA = 0; countA < allelePerLocus; countA++) { arrayGenotype[countN][countA] = vectAllelePool[countAPool]; countAPool++; } } // stat of year0 population // call statDiploidSingleLocus function statDiploidSingleLocus(arrayGenotype, allelePerLocus, nPopSize[0], nAllele, nOA, nEA, nHo, nHe, nF); // write to raw data array arrayRawOA[0][countI] = nOA; arrayRawEA[0][countI] = nEA; arrayRawHo[0][countI] = nHo; arrayRawHe[0][countI] = nHe; arrayRawF[0][countI] = nF; // fill age array for (countN = 0; countN < nPopSize[0]; countN++) { // age // when degreeOverlap = 0, if-test will always return true -> all individual start at age 0 // when degreeOverlap = 100, if-test will always return fales -> all individual start with random age value if (degreeOverlap < getRandomNumber(1, 100)) { arrayAge[countN] = 0; } else { arrayAge[countN] = getRandomNumber(0, (longevity - 1)); } } // start of year loop // each year loop has 2 major steps // 1. from previous year end to current year begin // generate survived population from the previous year // 2. from current year begin to current year end // generate breeding population of the current year // check age/replace old individuals reached longevity limit // geneerate new individuals if population growth occurred // // individual death events (old genotypes eliminated or replaced by new ones) can happen at both steps // death in step 1 is due to population decline (random death) // death in step 2 is due to reaching longevity limit (check age) // individual birth events always happen at step 2 // possible birth event1: old individuals reached longevity limit // possible birth event2: population growth occurred for (countY = 0; countY < nYear; countY++) { // check if asexual/monoecious or dioecious // if asexual/monoecious --> treat the population as a whole // else (dioecious) --> treat female and male separately if (reproSys >= 1 && reproSys <= 4) // asexual/monoecious, treat the population aas a whole { // step1: from previous year end to curent year begin // generate survived population (arrayGenotypeTemp) // check if population decline/stasis/growth // if decline -> get survivor list, copy survivor genotype and age // if stasis/growth -> copy all genotype and age from previous year end // check if population decline if (nPopSize[(countY+1)] < nPopSize[countY]) { nSurvivedTotal = nPopSize[(countY + 1)]; // fill population list vector for selecting surviving individuals for (countN = 0; countN < nPopSize[countY]; countN++) vectPopList.push_back(countN); // random suffle to generate survivor list random_shuffle(vectPopList.begin(), vectPopList.end()); // copy survivor genotype from arrayGenotype (previous year end) to arrayGenotypeTemp (current year begin) // copy survivor age (and increase by 1) for (countN = 0; countN < nSurvivedTotal; countN++) { arrayGenotypeTemp[countN][0] = arrayGenotype[(vectPopList[countN])][0]; arrayGenotypeTemp[countN][1] = arrayGenotype[(vectPopList[countN])][1]; arrayAgeTemp[countN] = (arrayAge[(vectPopList[countN])] + 1); } // end of copy genotype/age for loop // erase population list vector for next year vectPopList.erase(vectPopList.begin(), vectPopList.end()); } // end of if population decline else // population stasis or growth { nSurvivedTotal = nPopSize[countY]; // copy all genotype and age from previous year for (countN = 0; countN < nSurvivedTotal; countN++) { arrayGenotypeTemp[countN][0] = arrayGenotype[countN][0]; arrayGenotypeTemp[countN][1] = arrayGenotype[countN][1]; arrayAgeTemp[countN] = (arrayAge[countN] + 1); } } // end of population stasis or growth else // step2: from current year begin to current year end // generate breeding population // check individual age // if age reached longevity limit -> generate new individuals // check if population growth -> generate new individuals // copy genotype from arrayGenotypeTemp (current year begin) to arrayGenotype (current year end) // copy age from arrayAgeTemp to arrayAge // if population decline /stasis // -> total nPopSize[(countY+1)] individuals, all are from the previous year // if population growth // -> total nPopSize[(countY+1)] individuals, nPopSize[countY] are from the previous year // -> individuals from (nPopSize[countY]+1) to nPopSize[(countY+1)] are newborns of current year // generate breeding population // check if reached reproductive maturation, copy genotype to arrayParent nAdultTotal = 0; // reset number of adults for (countN = 0; countN < nSurvivedTotal; countN++) { if (arrayAgeTemp[countN] >= ageMaturation) { arrayParent[nAdultTotal][0] = arrayGenotypeTemp[countN][0]; arrayParent[nAdultTotal][1] = arrayGenotypeTemp[countN][1]; nAdultTotal++; } } // check survivor age // replace with new genotype if reached limit // else copy from arrayGenotypeTemp to arrayGenotype for (countN = 0; countN < nSurvivedTotal; countN++) { // check if reached longevity limit, generate new individual if (arrayAgeTemp[countN] >= longevity) { // reset age arrayAge[countN] = 0; // generate new genotype // total nAdultTotal breeders getNewDiSLGenotype(arrayParent, arrayNewGenotype, reproSys, nAdultTotal, nAdultFemale, parent1, parent2); // copy genotype arrayGenotype[countN][0] = arrayNewGenotype[0]; arrayGenotype[countN][1] = arrayNewGenotype[1]; } // end of check age if else { // copy age arrayAge[countN] = arrayAgeTemp[countN]; // copy genotype (from current year begin to current year end) arrayGenotype[countN][0] = arrayGenotypeTemp[countN][0]; arrayGenotype[countN][1] = arrayGenotypeTemp[countN][1]; } } // end of check age for loop // if population growth occurs // individual from (nPopSize[countY]) to (nPopSize[(countY+1)] - 1) are newborns of current year for (countN = nPopSize[countY]; countN < nPopSize[(countY+1)]; countN++) { // reset age (newborn age = 0) arrayAge[countN] = 0; // generate new genotype // total nAdultTotal breeders getNewDiSLGenotype(arrayParent, arrayNewGenotype, reproSys, nAdultTotal, nAdultFemale, parent1, parent2); // copy genotype arrayGenotype[countN][0] = arrayNewGenotype[0]; arrayGenotype[countN][1] = arrayNewGenotype[1]; } // end of population growth for loop } // end of asexual/monoecious if else // dioecious reproSys, treat female/male separately { // step1: from previous year end to curent year begin // generate breeding population (arrayGenotypeTemp) // check if nFemale/nMale decline/stasis/growth // if decline -> get survivor list, copy survivor genotype and age // if stasis/growth -> copy all genotype and age from previous year end // check if nFemale decline if (nFemale[(countY+1)] < nFemale[countY]) { // number of survived female = nFemale[(countY + 1)] nSurvivedFemale = nFemale[(countY + 1)]; // fill nFemale list vector for selecting surviving individuals for (countN = 0; countN < nFemale[countY]; countN++) vectPopList.push_back(countN); // random suffle to generate survivor list random_shuffle(vectPopList.begin(), vectPopList.end()); // copy survivor genotype from arrayGenotype (previous year end) to arrayGenotypeTemp (current year begin) // copy survivor age (and increase by 1) for (countN = 0; countN < nSurvivedFemale; countN++) { arrayGenotypeTemp[countN][0] = arrayGenotype[(vectPopList[countN])][0]; arrayGenotypeTemp[countN][1] = arrayGenotype[(vectPopList[countN])][1]; arrayAgeTemp[countN] = (arrayAge[(vectPopList[countN])] + 1); } // end of copy genotype/age for loop // erase population list vector for next year vectPopList.erase(vectPopList.begin(), vectPopList.end()); } // end of if nFemale decline else // nFemale stasis or growth { // number of survived female = nFemale[countY] nSurvivedFemale = nFemale[countY]; // copy all genotype and age from previous year for (countN = 0; countN < nSurvivedFemale; countN++) { arrayGenotypeTemp[countN][0] = arrayGenotype[countN][0]; arrayGenotypeTemp[countN][1] = arrayGenotype[countN][1]; arrayAgeTemp[countN] = (arrayAge[countN] + 1); } } // end of nFemale stasis or growth else // check if nMale decline if (nMale[(countY+1)] < nMale[countY]) { // number of survived male = nMale[(countY + 1)] nSurvivedMale = nMale[(countY + 1)]; // fill nMale list vector for selecting surviving individuals for (countN = 0; countN < nMale[countY]; countN++) vectPopList.push_back((countN + nFemale[countY])); // random suffle to generate survivor list random_shuffle(vectPopList.begin(), vectPopList.end()); // copy survivor genotype from arrayGenotype (previous year end) to arrayGenotypeTemp (current year begin) // copy survivor age (and increase by 1) for (countN = 0; countN < nSurvivedMale; countN++) { arrayGenotypeTemp[(countN + nSurvivedFemale)][0] = arrayGenotype[(vectPopList[countN])][0]; arrayGenotypeTemp[(countN + nSurvivedFemale)][1] = arrayGenotype[(vectPopList[countN])][1]; arrayAgeTemp[(countN + nSurvivedFemale)] = (arrayAge[(vectPopList[countN])] + 1); } // end of copy genotype/age for loop // erase population list vector for next year vectPopList.erase(vectPopList.begin(), vectPopList.end()); } // end of if nMale decline else // nMale stasis or growth { // number of survived male = nMale[countY] nSurvivedMale = nMale[countY]; // copy all genotype and age from previous year for (countN = 0; countN < nSurvivedMale; countN++) { arrayGenotypeTemp[(countN + nSurvivedFemale)][0] = arrayGenotype[(countN + nFemale[countY])][0]; arrayGenotypeTemp[(countN + nSurvivedFemale)][1] = arrayGenotype[(countN + nFemale[countY])][1]; arrayAgeTemp[(countN + nSurvivedFemale)] = (arrayAge[(countN + nFemale[countY])] + 1); } } // end of nMale stasis or growth else // step2: from current year begin to current year end // generate breeding population // check individual age // if age reached longevity limit -> generate new individuals // check if population growth -> generate new individuals // copy genotype from arrayGenotypeTemp (current year begin) to arrayGenotype (current year end) // copy age from arrayAgeTemp to arrayAge // generate female/male newborns if needed // generate breeding population nAdultFemale = 0; for (countN = 0; countN < nSurvivedFemale; countN++) { if (arrayAgeTemp[countN] >= ageMaturation) { arrayParent[nAdultFemale][0] = arrayGenotypeTemp[countN][0]; arrayParent[nAdultFemale][1] = arrayGenotypeTemp[countN][1]; nAdultFemale++; } } nAdultMale = 0; for (countN = 0; countN < nSurvivedMale; countN++) { if (arrayAgeTemp[(countN + nSurvivedFemale)] >= ageMaturation) { arrayParent[(nAdultMale + nAdultFemale)][0] = arrayGenotypeTemp[(countN + nSurvivedFemale)][0]; arrayParent[(nAdultMale + nAdultFemale)][1] = arrayGenotypeTemp[(countN + nSurvivedFemale)][1]; nAdultMale++; } } nAdultTotal = (nAdultFemale + nAdultMale); // get parent id for diecious model with single reproducing male/pair // single reproducing male // get male id if (reproSys == 6) { parent2 = getRandomNumber(nAdultFemale, (nAdultTotal - 1)); } // single reproducing pair // get both parent ids if (reproSys == 7) { parent1 = getRandomNumber(0, (nAdultFemale - 1)); parent2 = getRandomNumber(nAdultFemale, (nAdultTotal - 1)); } // check survived female age // replace with new genotype if reached limit // else copy from arrayGenotypeTemp to arrayGenotype for (countN = 0; countN < nSurvivedFemale; countN++) { // check if reached longevity limit, generate new individual if (arrayAgeTemp[countN] >= longevity) { // reset age arrayAge[countN] = 0; // generate new genotype // total nAdultTotal breeders getNewDiSLGenotype(arrayParent, arrayNewGenotype, reproSys, nAdultTotal, nAdultFemale, parent1, parent2); // copy genotype arrayGenotype[countN][0] = arrayNewGenotype[0]; arrayGenotype[countN][1] = arrayNewGenotype[1]; } // end of check age if else { // copy age arrayAge[countN] = arrayAgeTemp[countN]; // copy genotype (from current year begin to current year end) arrayGenotype[countN][0] = arrayGenotypeTemp[countN][0]; arrayGenotype[countN][1] = arrayGenotypeTemp[countN][1]; } } // end of check age for loop // if nFemale growth occurs // individual from (nFemale[countY]) to (nFemale[(countY+1)] - 1) are newborns of current year for (countN = nFemale[countY]; countN < nFemale[(countY+1)]; countN++) { // reset age (newborn age = 0) arrayAge[countN] = 0; // generate new genotype // total nAdultTotal breeders getNewDiSLGenotype(arrayParent, arrayNewGenotype, reproSys, nAdultTotal, nAdultFemale, parent1, parent2); // copy genotype arrayGenotype[countN][0] = arrayNewGenotype[0]; arrayGenotype[countN][1] = arrayNewGenotype[1]; } // end of population growth for loop // check survived male age // replace with new genotype if reached limit // else copy from arrayGenotypeTemp to arrayGenotype for (countN = 0; countN < nSurvivedMale; countN++) { // check if reached longevity limit, generate new individual if (arrayAgeTemp[(countN + nSurvivedFemale)] >= longevity) { // reset age arrayAge[(countN + nFemale[(countY + 1)])] = 0; // generate new genotype // total nAdultTotal breeders getNewDiSLGenotype(arrayParent, arrayNewGenotype, reproSys, nAdultTotal, nAdultFemale, parent1, parent2); // copy genotype arrayGenotype[(countN + nFemale[(countY + 1)])][0] = arrayNewGenotype[0]; arrayGenotype[(countN + nFemale[(countY + 1)])][1] = arrayNewGenotype[1]; } // end of check age if else { // copy age arrayAge[(countN + nFemale[(countY + 1)])] = arrayAgeTemp[(countN + nSurvivedFemale)]; // copy genotype (from current year begin to current year end) arrayGenotype[(countN + nFemale[(countY + 1)])][0] = arrayGenotypeTemp[(countN + nSurvivedFemale)][0]; arrayGenotype[(countN + nFemale[(countY + 1)])][1] = arrayGenotypeTemp[(countN + nSurvivedFemale)][1]; } } // end of check age for loop // if nMale growth occurs for (countN = 0; countN < (nMale[(countY+1)] - nMale[countY]); countN++) { // reset age (newborn age = 0) arrayAge[(countN + nFemale[(countY + 1)] + nSurvivedMale)] = 0; // generate new genotype // total nAdultTotal breeders getNewDiSLGenotype(arrayParent, arrayNewGenotype, reproSys, nAdultTotal, nAdultFemale, parent1, parent2); // copy genotype arrayGenotype[(countN + nFemale[(countY + 1)] + nSurvivedMale)][0] = arrayNewGenotype[0]; arrayGenotype[(countN + nFemale[(countY + 1)] + nSurvivedMale)][1] = arrayNewGenotype[1]; } // end of nMale growth for loop } // end of dioecious reproSys else // stat // call statDiploidSingleLocus function statDiploidSingleLocus(arrayGenotype, allelePerLocus, nPopSize[(countY+1)], nAllele, nOA, nEA, nHo, nHe, nF); // write to data array arrayRawOA[(countY + 1)][countI] = nOA; arrayRawEA[(countY + 1)][countI] = nEA; arrayRawHo[(countY + 1)][countI] = nHo; arrayRawHe[(countY + 1)][countI] = nHe; arrayRawF[(countY + 1)][countI] = nF; } // end of year loop // check boolOutputRawGenotype if (boolOutputRawGenotype == 1) { // call writeDiSLRawGenotype function writeDiSLRawGenotype(rawGenotypeOutputFile, arrayGenotype, allelePerLocus, nPopSize[nYear]); } cout << '.'; // use dot to show iteration progress } // end of iteration loop // display message cout << endl; cout << "Simulation completed." << endl; cout << endl; cout << "Preparing output file..." << endl; cout << endl; // write simulation settings to output file // call writeVPSSimSetting function writeVPSSimSetting(outputFile, degreeOverlap, reproSys, longevity, ageMaturation, nYear, nIteration); // write summary output file writeDiSLSummary(outputFile, nYear, nIteration, arrayRawOA, arrayRawEA, arrayRawHo, arrayRawHe, arrayRawF); // call closeFile function // close input file closeInputFile(alleleFreqInput); closeInputFile(popSizeInput); // close output file closeOutputFile(outputFile); // check boolOutputRawGenotype if (boolOutputRawGenotype == 1) { closeOutputFile(rawGenotypeOutputFile); } } // end of simDiploidSingleLocusVPS function // simDiploidMultiLocusCPS function // perform diploid multilocus constant population size simulation void simDiploidMultiLocusCPS(void) { // constants const int allelePerLocus = 2; // number of alleles per locus // declare loop counting variables int countN = 0; // count individual int countL = 0; // count locus int countA = 0; // count allele int countAId = 0; // count allele id int countY = 0; // count year int countI = 0; // count iteration // file operation // input file - allele freq setting // declare fileName fileName genotypeInputName; // input - genotype data // declare input file ifstream genotypeInputFile; // call openInput function openInput(genotypeInputFile, genotypeInputName); // output file - simulation summary // declare fileName fileName outputFileName; // output - simulation summary // declare output file ofstream outputFile; // call openOutput function openOutput(outputFile, outputFileName); // call writeOutputHeader function // write header to output file writeOutputHeader(outputFile); // output file - raw genotype data // declare fileName fileName rawGenotypeOutputName; // declare raw data output file ofstream rawGenotypeOutputFile; // read input file // declare variables int nSize = 0; // sample size (# of individuals) int nLoci = 0; // # of loci // display cout << endl; cout << "Reading input file..." << endl; // read nSize nSize = readInputFileNSize(genotypeInputFile); // read nLoci nLoci = readInputFileNLoci(genotypeInputFile); // initiate input array // array to hold locus name locusName *arrayInputLocus; arrayInputLocus = new locusName[nLoci]; if (arrayInputLocus == NULL) memoryError(); // array to hold individual ID individualId *arrayInputId; arrayInputId = new individualId[nSize]; if (arrayInputId == NULL) memoryError(); // array to hold input genotype data char ***arrayInput; arrayInput = new char** [nSize]; if (arrayInput == NULL) memoryError(); for (countN = 0; countN < nSize; countN++) { arrayInput[countN] = new char* [nLoci]; if (arrayInput == NULL) memoryError(); for (countL = 0; countL < nLoci; countL++) { arrayInput[countN][countL] = new char [allelePerLocus]; if (arrayInput == NULL) memoryError(); } } // read input data // locus name readInputFileLocusName(genotypeInputFile, arrayInputLocus, nLoci); // genotype readInputFileGenotype(genotypeInputFile, arrayInputId, arrayInput, nSize, nLoci, allelePerLocus); // display imported input data displayGenotypeInput(arrayInputLocus, arrayInputId, arrayInput, nSize, nLoci, allelePerLocus); // calculate allele freq from input file // initiate data array // arrayAlleleCount // array for allele count [nLoci][maxAllele] int **arrayAlleleCount; arrayAlleleCount = new int* [nLoci] ; if (arrayAlleleCount == NULL) memoryError(); for (countL = 0; countL < nLoci; countL++) { arrayAlleleCount[countL] = new int [maxAllele]; if (arrayAlleleCount == NULL) memoryError(); } // alleleTotalCount // array for allele total count [nLoci] int *alleleTotalCount; alleleTotalCount = new int [nLoci]; if (alleleTotalCount == NULL) memoryError(); // arrayAlleleFreq // array for allele freq [nLoci][maxAllele] float **arrayAlleleFreq; arrayAlleleFreq = new float* [nLoci] ; if (arrayAlleleFreq == NULL) memoryError(); for (countL = 0; countL < nLoci; countL++) { arrayAlleleFreq[countL] = new float [maxAllele]; if (arrayAlleleFreq == NULL) memoryError(); } // call statGenotypeInput function statGenotypeInput(arrayInput, arrayAlleleCount, alleleTotalCount, arrayAlleleFreq, nSize, nLoci, allelePerLocus); // write input allele freq data to output file // call writeAlleleFreq function writeAlleleFreq(outputFile, nLoci, arrayInputLocus, arrayAlleleFreq); // get simulation settings // declare variables int degreeOverlap; // generation overlapping setting int reproSys; // reproductive system option int longevity; // longevity of the organism int ageMaturation; // reproductive maturation age of the organism int nPB; // pre-bottleneck population size int nB; // bottlneck population size // number of female/male for diecious models int nFemale; int nMale; int nYear; // duration of bottleneck int nIteration; // number of iterations of the run int boolOutputRawGenotype; // bool variable for output raw genotype data // call getCPSSimSetting function getCPSSimSetting(degreeOverlap, reproSys, longevity, ageMaturation, nPB, nB, nFemale, nMale, nYear, nIteration, boolOutputRawGenotype); // check boolOutputRawGenotype if (boolOutputRawGenotype == 1) { // open raw data output file // call openRawGenotypeOutput function openRawGenotypeOutput(rawGenotypeOutputFile, rawGenotypeOutputName); // write header writeMLRawGenotypeHeader(rawGenotypeOutputFile, arrayInputLocus, nLoci); } // declare simulation variables int parent1 = 0; // id of parent1 int parent2 = 0; // id of parent2 int nSurvivedFemale = 0; // number of survived females int nSurvivedMale = 0; // number of survived males int nSurvivedTotal = 0; // number of survivors int nAdultFemale = 0; // number of adult females int nAdultMale = 0; // number of adult Males int nAdultTotal = 0; // number of adults // initiate data arrays // pre-bottleneck population genotype [nPB][nLoci][allelePerLocus] int ***arrayPBP; arrayPBP = new int** [nPB]; if (arrayPBP == NULL) memoryError(); for (countN = 0; countN < nPB; countN++) { arrayPBP[countN] = new int* [nLoci]; if (arrayPBP == NULL) memoryError(); for (countL = 0; countL < nLoci; countL++) { arrayPBP[countN][countL] = new int[allelePerLocus]; if (arrayPBP == NULL) memoryError(); } } // bottleneck population genotype [nB][nLoci][allelePerLocus] int ***arrayBP; arrayBP = new int** [nB]; if (arrayBP == NULL) memoryError(); for (countN = 0; countN < nB; countN++) { arrayBP[countN] = new int* [nLoci]; if (arrayBP == NULL) memoryError(); for (countL = 0; countL < nLoci; countL++) { arrayBP[countN][countL] = new int[allelePerLocus]; if (arrayBP == NULL) memoryError(); } } // breeding population genotype [nB][nLoci][allelePerLocus] int ***arrayParent; arrayParent = new int** [nB]; if (arrayParent == NULL) memoryError(); for (countN = 0; countN < nB; countN++) { arrayParent[countN] = new int* [nLoci]; if (arrayParent == NULL) memoryError(); for (countL = 0; countL < nLoci; countL++) { arrayParent[countN][countL] = new int[allelePerLocus]; if (arrayParent == NULL) memoryError(); } } // genotype array for new individuals int **arrayNewGenotype; arrayNewGenotype = new int* [nLoci]; if (arrayNewGenotype == NULL) memoryError(); for (countL = 0; countL < nLoci; countL++) { arrayNewGenotype[countL] = new int [allelePerLocus]; if (arrayNewGenotype == NULL) memoryError(); } // age [nB] int *arrayAge; arrayAge = new int [nB]; if (arrayAge == NULL) memoryError(); // initiate stat array [nLoci] // OA (observed number of alleles) int *arrayOA; arrayOA = new int [nLoci]; if (arrayOA == NULL ) memoryError(); // EA (effective number of alleles) float *arrayEA; arrayEA = new float [nLoci]; if (arrayEA == NULL) memoryError(); // Ho (observed heterozygosity) float *arrayHo; arrayHo = new float [nLoci]; if (arrayHo == NULL) memoryError(); // He (expected heterozygosity) float *arrayHe; arrayHe = new float [nLoci]; if (arrayHe == NULL) memoryError(); // F (fixation index) float *arrayF; arrayF = new float [nLoci]; if (arrayF == NULL) memoryError(); // initiate raw data array [nLoci][nYear][nIteration] // OA [nLoci][nYear][nIteration] int ***arrayRawOA; arrayRawOA = new int** [nLoci]; if (arrayRawOA == NULL) memoryError(); for (countL = 0; countL < nLoci; countL++) { arrayRawOA[countL] = new int* [(nYear + 1)]; if (arrayRawOA == NULL) memoryError(); for (countY = 0; countY < (nYear + 1); countY++) { arrayRawOA[countL][countY] = new int [nIteration]; if (arrayRawOA == NULL) memoryError(); } } // EA [nLoci][nYear][nIteration] float ***arrayRawEA; arrayRawEA = new float** [nLoci]; if (arrayRawEA == NULL) memoryError(); for (countL = 0; countL < nLoci; countL++) { arrayRawEA[countL] = new float* [(nYear + 1)]; if (arrayRawEA == NULL) memoryError(); for (countY = 0; countY < (nYear + 1); countY++) { arrayRawEA[countL][countY] = new float [nIteration]; if (arrayRawEA == NULL) memoryError(); } } // Ho [nLoci][nYear][nIteration] float ***arrayRawHo; arrayRawHo = new float** [nLoci]; if (arrayRawHo == NULL) memoryError(); for (countL = 0; countL < nLoci; countL++) { arrayRawHo[countL] = new float* [(nYear + 1)]; if (arrayRawHo == NULL) memoryError(); for (countY = 0; countY < (nYear + 1); countY++) { arrayRawHo[countL][countY] = new float [nIteration]; if (arrayRawHo == NULL) memoryError(); } } // He [nLoci][nYear][nIteration] float ***arrayRawHe; arrayRawHe = new float** [nLoci]; if (arrayRawHe == NULL) memoryError(); for (countL = 0; countL < nLoci; countL++) { arrayRawHe[countL] = new float* [(nYear + 1)]; if (arrayRawHe == NULL) memoryError(); for (countY = 0; countY < (nYear + 1); countY++) { arrayRawHe[countL][countY] = new float [nIteration]; if (arrayRawHe == NULL) memoryError(); } } // F [nLoci][nYear][nIteration] float ***arrayRawF; arrayRawF = new float** [nLoci]; if (arrayRawF == NULL) memoryError(); for (countL = 0; countL < nLoci; countL++) { arrayRawF[countL] = new float* [(nYear + 1)]; if (arrayRawF == NULL) memoryError(); for (countY = 0; countY < (nYear + 1); countY++) { arrayRawF[countL][countY] = new float [nIteration]; if (arrayRawF == NULL) memoryError(); } } // vector for list // allele pool for generating the pre-bottleneck population vector vectAllelePool; // pre bottleneck population list [nPB] // list of all individuals in the pre bottleneck population // for selecting individuals survived bottleneck vector vectPBPList; for (countN = 0; countN < nPB; countN++) vectPBPList.push_back(countN); // start simulation // start of iteration loop for (countI = 0; countI < nIteration; countI++) { // generate the pre bottleneck population for (countL = 0; countL < nLoci; countL++) { // reset allele pool vector vectAllelePool.erase(vectAllelePool.begin(), vectAllelePool.end()); // fill allele pool vector with the number of allele counts from input data for (countA = 0; countA < maxAllele; countA++) { for (countAId = 0; countAId < arrayAlleleCount[countL][countA]; countAId++) { vectAllelePool.push_back(countA); } } // fill pre bottleneck population genotype array for (countN = 0; countN < nPB; countN++) { for (countA = 0; countA < allelePerLocus; countA++) { // random shuffle allele pool random_shuffle(vectAllelePool.begin(), vectAllelePool.end()); // get allele from allele pool arrayPBP[countN][countL][countA] = vectAllelePool[0]; } } } // end of countL loop for generation pre bottleneck population // stat of pre bottleneck population // call statDiploidMultiLocus function statDiploidMultiLocus(arrayPBP, allelePerLocus, nPB, nLoci, arrayOA, arrayEA, arrayHo, arrayHe, arrayF); // write to raw data array for (countL = 0; countL < nLoci; countL++) { arrayRawOA[countL][0][countI] = arrayOA[countL]; arrayRawEA[countL][0][countI] = arrayEA[countL]; arrayRawHo[countL][0][countI] = arrayHo[countL]; arrayRawHe[countL][0][countI] = arrayHe[countL]; arrayRawF[countL][0][countI] = arrayF[countL]; } // generate bottleneck population random_shuffle(vectPBPList.begin(), vectPBPList.end()); // survivor list for (countN = 0; countN < nB; countN++) { // genotype for (countL = 0; countL < nLoci; countL++) { for (countA = 0; countA < allelePerLocus; countA++) { arrayBP[countN][countL][countA] = arrayPBP[(vectPBPList[countN])][countL][countA]; } } // age // when degreeOverlap = 0, if-test will always return true -> all individual start at age 0 // when degreeOverlap = 100, if-test will always return fales -> all individual start with random age value if (degreeOverlap < getRandomNumber(1, 100)) { arrayAge[countN] = 0; } else { arrayAge[countN] = getRandomNumber(0, (longevity - 1)); } } // start of year loop // each year loop has 2 major steps // 1. from previous year end to current year begin // generate survived population from the previous year // 2. from current year begin to current year end // generate breeding population of the current year // check age/replace old individuals reached longevity limit // geneerate new individuals if population growth occurred // // individual death events (old genotypes eliminated or replaced by new ones) can happen at both steps // death in step 1 is due to population decline (random death) // death in step 2 is due to reaching longevity limit (check age) // individual birth events always happen at step 2 // possible birth event1: old individuals reached longevity limit // possible birth event2: population growth occurred for (countY = 0; countY < nYear; countY++) { // check if asexual/monoecious or dioecious // if asexual/monoecious --> treat the population as a whole // else (dioecious) --> treat female and male separately if (reproSys >= 1 && reproSys <= 4) // asexual/monoecious, treat the population aas a whole { // step1: from previous year end to curent year begin // in constant population size module all individuals survived // do nothing to genotype array nSurvivedTotal = nB; // age increased by 1 for (countN = 0; countN < nSurvivedTotal; countN++) { arrayAge[countN]++; } // step2: from current year begin to current year end // generate breeding population // check individual age // if age reached longevity limit -> generate new individuals // generate breeding population // check if reached reproductive maturation, copy genotype to arrayParent nAdultTotal = 0; // reset number of adults for (countN = 0; countN < nSurvivedTotal; countN++) { if (arrayAge[countN] >= ageMaturation) { for (countL = 0; countL < nLoci; countL++) { arrayParent[nAdultTotal][countL][0] = arrayBP[countN][countL][0]; arrayParent[nAdultTotal][countL][1] = arrayBP[countN][countL][1]; } nAdultTotal++; } } // check survivor age // reset age and replace with new genotype if reached limit // else do nothing to genotype and age array for (countN = 0; countN < nSurvivedTotal; countN++) { // check if reached longevity limit, generate new individual if (arrayAge[countN] >= longevity) { // reset age arrayAge[countN] = 0; // generate new genotype // total nAdultTotal breeders getNewDiMLGenotype(arrayParent, arrayNewGenotype, reproSys, nAdultTotal, nAdultFemale, parent1, parent2, nLoci); // copy genotype for (countL = 0; countL < nLoci; countL++) { arrayBP[countN][countL][0] = arrayNewGenotype[countL][0]; arrayBP[countN][countL][1] = arrayNewGenotype[countL][1]; } } // end of check age if // else do nothing to genotype and age array } // end of check age for loop } // end of asexual/monoecious if else // dioecious reproSys, treat female/male separately { // step1: from previous year end to curent year begin // in constant population size module all individuals survived // do nothing to genotype array nSurvivedFemale = nFemale; nSurvivedMale = nMale; nSurvivedTotal = nB; // age increased by 1 for (countN = 0; countN < nSurvivedTotal; countN++) { arrayAge[countN]++; } // step2: from current year begin to current year end // generate breeding population // check individual age // if age reached longevity limit -> generate new individuals // generate breeding population // check if reached reproductive maturation, copy genotype to arrayParent // generate breeding population nAdultFemale = 0; for (countN = 0; countN < nSurvivedFemale; countN++) { if (arrayAge[countN] >= ageMaturation) { for (countL = 0; countL < nLoci; countL++) { arrayParent[nAdultFemale][countL][0] = arrayBP[countN][countL][0]; arrayParent[nAdultFemale][countL][1] = arrayBP[countN][countL][1]; } nAdultFemale++; } } nAdultMale = 0; for (countN = 0; countN < nSurvivedMale; countN++) { if (arrayAge[(countN + nSurvivedFemale)] >= ageMaturation) { for (countL = 0; countL < nLoci; countL++) { arrayParent[(nAdultMale + nAdultFemale)][countL][0] = arrayBP[(countN + nSurvivedFemale)][countL][0]; arrayParent[(nAdultMale + nAdultFemale)][countL][1] = arrayBP[(countN + nSurvivedFemale)][countL][1]; } nAdultMale++; } } nAdultTotal = (nAdultFemale + nAdultMale); // get parent id for diecious model with single reproducing male/pair // single reproducing male // get male id if (reproSys == 6) { parent2 = getRandomNumber(nAdultFemale, (nAdultTotal - 1)); } // single reproducing pair // get both parent ids if (reproSys == 7) { parent1 = getRandomNumber(0, (nAdultFemale - 1)); parent2 = getRandomNumber(nAdultFemale, (nAdultTotal - 1)); } // check survivor age // reset age and replace with new genotype if reached limit // else do nothing to genotype and age array for (countN = 0; countN < nSurvivedTotal; countN++) { // check if reached longevity limit, generate new individual if (arrayAge[countN] >= longevity) { // reset age arrayAge[countN] = 0; // generate new genotype // total nAdultTotal breeders getNewDiMLGenotype(arrayParent, arrayNewGenotype, reproSys, nAdultTotal, nAdultFemale, parent1, parent2, nLoci); // copy genotype for (countL = 0; countL < nLoci; countL++) { arrayBP[countN][countL][0] = arrayNewGenotype[countL][0]; arrayBP[countN][countL][1] = arrayNewGenotype[countL][1]; } } // end of check age if // else do nothing to genotype and age array } // end of check age for loop } // end of dioecious repro sys else // stat of bottleneck population // call statDiploidMultiLocus function statDiploidMultiLocus(arrayBP, allelePerLocus, nB, nLoci, arrayOA, arrayEA, arrayHo, arrayHe, arrayF); // write to raw data array for (countL = 0; countL < nLoci; countL++) { arrayRawOA[countL][(countY + 1)][countI] = arrayOA[countL]; arrayRawEA[countL][(countY + 1)][countI] = arrayEA[countL]; arrayRawHo[countL][(countY + 1)][countI] = arrayHo[countL]; arrayRawHe[countL][(countY + 1)][countI] = arrayHe[countL]; arrayRawF[countL][(countY + 1)][countI] = arrayF[countL]; } } // end of year loop // check boolOutputRawGenotype if (boolOutputRawGenotype == 1) { // call writeDiSLRawGenotype function writeDiMLRawGenotype(rawGenotypeOutputFile, arrayBP, allelePerLocus, nB, nLoci); } cout << '.'; // use dot to show iteration progress } // end of iteration loop // display message cout << endl; cout << "Simulation completed." << endl; cout << endl; cout << "Preparing output file..." << endl; cout << endl; // write simulation settings to output file // call writeCPSSimSetting function writeCPSSimSetting(outputFile, degreeOverlap, reproSys, longevity, ageMaturation, nPB, nB, nFemale, nMale, nYear, nIteration); // write summary output file writeDiMLSummary(outputFile, nLoci, arrayInputLocus, nYear, nIteration, arrayRawOA, arrayRawEA, arrayRawHo, arrayRawHe, arrayRawF); // call closeFile function // close input file closeInputFile(genotypeInputFile); // close output file closeOutputFile(outputFile); // check boolOutputRawGenotype if (boolOutputRawGenotype == 1) { closeOutputFile(rawGenotypeOutputFile); } } // end of diploidMultiLocusCPS function // simDiploidMultiLocusVPS function // perform diploid multilocus variable population size simulation void simDiploidMultiLocusVPS(void) { // constants const int allelePerLocus = 2; // number of alleles per locus // declare loop counting variables int countN = 0; // count individual int countL = 0; // count locus int countA = 0; // count allele int countAId = 0; // count allele id int countY = 0; // count year int countI = 0; // count iteration // file operation // input file - allele freq setting // declare fileName fileName genotypeInputName; // input - genotype data // declare input file ifstream genotypeInputFile; // call openInput function openInput(genotypeInputFile, genotypeInputName); // input file - population size setting // declare fileName fileName popSizeInputName; // input - population size // declare input file ifstream popSizeInput; // output file - simulation summary // declare fileName fileName outputFileName; // output - simulation summary // declare output file ofstream outputFile; // call openOutput function openOutput(outputFile, outputFileName); // call writeOutputHeader function // write header to output file writeOutputHeader(outputFile); // output file - raw genotype data // declare fileName fileName rawGenotypeOutputName; // declare raw data output file ofstream rawGenotypeOutputFile; // read input file // declare variables int nSize = 0; // sample size (# of individuals) int nLoci = 0; // # of loci // display cout << endl; cout << "Reading input file..." << endl; // read nSize nSize = readInputFileNSize(genotypeInputFile); // read nLoci nLoci = readInputFileNLoci(genotypeInputFile); // initiate input array // array to hold locus name locusName *arrayInputLocus; arrayInputLocus = new locusName[nLoci]; if (arrayInputLocus == NULL) memoryError(); // array to hold individual ID individualId *arrayInputId; arrayInputId = new individualId[nSize]; if (arrayInputId == NULL) memoryError(); // array to hold input genotype data char ***arrayInput; arrayInput = new char** [nSize]; if (arrayInput == NULL) memoryError(); for (countN = 0; countN < nSize; countN++) { arrayInput[countN] = new char* [nLoci]; if (arrayInput == NULL) memoryError(); for (countL = 0; countL < nLoci; countL++) { arrayInput[countN][countL] = new char [allelePerLocus]; if (arrayInput == NULL) memoryError(); } } // read input data // locus name readInputFileLocusName(genotypeInputFile, arrayInputLocus, nLoci); // genotype readInputFileGenotype(genotypeInputFile, arrayInputId, arrayInput, nSize, nLoci, allelePerLocus); // display imported input data displayGenotypeInput(arrayInputLocus, arrayInputId, arrayInput, nSize, nLoci, allelePerLocus); // calculate allele freq from input file // initiate data array // arrayAlleleCount // array for allele count [nLoci][maxAllele] int **arrayAlleleCount; arrayAlleleCount = new int* [nLoci] ; if (arrayAlleleCount == NULL) memoryError(); for (countL = 0; countL < nLoci; countL++) { arrayAlleleCount[countL] = new int [maxAllele]; if (arrayAlleleCount == NULL) memoryError(); } // alleleTotalCount // array for allele total count [nLoci] int *alleleTotalCount; alleleTotalCount = new int [nLoci]; if (alleleTotalCount == NULL) memoryError(); // arrayAlleleFreq // array for allele freq [nLoci][maxAllele] float **arrayAlleleFreq; arrayAlleleFreq = new float* [nLoci] ; if (arrayAlleleFreq == NULL) memoryError(); for (countL = 0; countL < nLoci; countL++) { arrayAlleleFreq[countL] = new float [maxAllele]; if (arrayAlleleFreq == NULL) memoryError(); } // call statGenotypeInput function statGenotypeInput(arrayInput, arrayAlleleCount, alleleTotalCount, arrayAlleleFreq, nSize, nLoci, allelePerLocus); // write input allele freq data to output file // call writeAlleleFreq function writeAlleleFreq(outputFile, nLoci, arrayInputLocus, arrayAlleleFreq); // get simulation settings // declare variables int degreeOverlap; // generation overlapping setting int reproSys; // reproductive system option int longevity; // longevity of the organism int ageMaturation; // reproductive maturation age of the organism int nIteration; // number of iterations of the run int boolOutputRawGenotype; // bool variable for output raw genotype data // call getVPSSimSetting function getVPSSimSetting(degreeOverlap, reproSys, longevity, ageMaturation, nIteration, boolOutputRawGenotype); // check boolOutputRawGenotype if (boolOutputRawGenotype == 1) { // open raw data output file // call openRawGenotypeOutput function openRawGenotypeOutput(rawGenotypeOutputFile, rawGenotypeOutputName); // write header writeMLRawGenotypeHeader(rawGenotypeOutputFile, arrayInputLocus, nLoci); } // get population size settings // call openPopSizeInput function openPopSizeInput(popSizeInput, popSizeInputName); // read population size setting input file // declare variables int nYear; // duration of bottleneck // display cout << endl; cout << "Reading population size setting input file..." << endl; // read nYear nYear = readInputFileNYear(popSizeInput); // initiate pop size input data array int *nPopSize; // array for popsize [(nYear+1)] nPopSize = new int [(nYear+1)]; if (nPopSize == NULL) memoryError(); int *nFemale; // array for nFemale [(nYear+1)] nFemale = new int [(nYear+1)]; if (nFemale == NULL) memoryError(); int *nMale; // array for nMale [(nYear+1)] nMale = new int [(nYear+1)]; if (nMale == NULL) memoryError(); // read population size input readInputFilePopSize(popSizeInput, nPopSize, nFemale, nMale, nYear, reproSys); // write imported population size input data writeVPopSizeInput(outputFile, nPopSize, nFemale, reproSys, nYear); // display imported population size input data displayVPopSizeInput(nPopSize, nFemale, reproSys, nYear); // check population size input checkVPopSsizeInput(nPopSize, nFemale, reproSys, nYear); // get max pop size int nMaxPopSize; // maximum population size in pop size input nMaxPopSize = getNMaxPopSize(nPopSize, nYear); // declare simulation variables int parent1 = 0; // id of parent1 int parent2 = 0; // id of parent2 int nSurvivedFemale = 0; // number of survived females int nSurvivedMale = 0; // number of survived males int nSurvivedTotal = 0; // number of survivors int nAdultFemale = 0; // number of adult females int nAdultMale = 0; // number of adult Males int nAdultTotal = 0; // number of adults // initiate data arrays // population genotype [nMaxPopSize][nLoci][allelePerLocus] int ***arrayGenotype; arrayGenotype = new int** [nMaxPopSize]; if (arrayGenotype == NULL) memoryError(); for (countN = 0; countN < nMaxPopSize; countN++) { arrayGenotype[countN] = new int* [nLoci]; if (arrayGenotype == NULL) memoryError(); for (countL = 0; countL < nLoci; countL++) { arrayGenotype[countN][countL] = new int[allelePerLocus]; if (arrayGenotype == NULL) memoryError(); } } // population genotype (temp) [nMaxPopSize][nLoci][allelePerLocus] int ***arrayGenotypeTemp; arrayGenotypeTemp = new int** [nMaxPopSize]; if (arrayGenotypeTemp == NULL) memoryError(); for (countN = 0; countN < nMaxPopSize; countN++) { arrayGenotypeTemp[countN] = new int* [nLoci]; if (arrayGenotypeTemp == NULL) memoryError(); for (countL = 0; countL < nLoci; countL++) { arrayGenotypeTemp[countN][countL] = new int[allelePerLocus]; if (arrayGenotypeTemp == NULL) memoryError(); } } // breeding population genotype [nMaxPopSize][nLoci][allelePerLocus] int ***arrayParent; arrayParent = new int** [nMaxPopSize]; if (arrayParent == NULL) memoryError(); for (countN = 0; countN < nMaxPopSize; countN++) { arrayParent[countN] = new int* [nLoci]; if (arrayParent == NULL) memoryError(); for (countL = 0; countL < nLoci; countL++) { arrayParent[countN][countL] = new int[allelePerLocus]; if (arrayParent == NULL) memoryError(); } } // genotype array for new individuals int **arrayNewGenotype; arrayNewGenotype = new int* [nLoci]; if (arrayNewGenotype == NULL) memoryError(); for (countL = 0; countL < nLoci; countL++) { arrayNewGenotype[countL] = new int [allelePerLocus]; if (arrayNewGenotype == NULL) memoryError(); } // age [nMaxPopSize] int *arrayAge; arrayAge = new int [nMaxPopSize]; if (arrayAge == NULL) memoryError(); // age (temp) [nMaxPopSize] int *arrayAgeTemp; arrayAgeTemp = new int [nMaxPopSize]; if (arrayAgeTemp == NULL) memoryError(); // initiate stat array [nLoci] // OA (observed number of alleles) int *arrayOA; arrayOA = new int [nLoci]; if (arrayOA == NULL ) memoryError(); // EA (effective number of alleles) float *arrayEA; arrayEA = new float [nLoci]; if (arrayEA == NULL) memoryError(); // Ho (observed heterozygosity) float *arrayHo; arrayHo = new float [nLoci]; if (arrayHo == NULL) memoryError(); // He (expected heterozygosity) float *arrayHe; arrayHe = new float [nLoci]; if (arrayHe == NULL) memoryError(); // F (fixation index) float *arrayF; arrayF = new float [nLoci]; if (arrayF == NULL) memoryError(); // initiate raw data array [nLoci][nYear][nIteration] // OA [nLoci][nYear][nIteration] int ***arrayRawOA; arrayRawOA = new int** [nLoci]; if (arrayRawOA == NULL) memoryError(); for (countL = 0; countL < nLoci; countL++) { arrayRawOA[countL] = new int* [(nYear + 1)]; if (arrayRawOA == NULL) memoryError(); for (countY = 0; countY < (nYear + 1); countY++) { arrayRawOA[countL][countY] = new int [nIteration]; if (arrayRawOA == NULL) memoryError(); } } // EA [nLoci][nYear][nIteration] float ***arrayRawEA; arrayRawEA = new float** [nLoci]; if (arrayRawEA == NULL) memoryError(); for (countL = 0; countL < nLoci; countL++) { arrayRawEA[countL] = new float* [(nYear + 1)]; if (arrayRawEA == NULL) memoryError(); for (countY = 0; countY < (nYear + 1); countY++) { arrayRawEA[countL][countY] = new float [nIteration]; if (arrayRawEA == NULL) memoryError(); } } // Ho [nLoci][nYear][nIteration] float ***arrayRawHo; arrayRawHo = new float** [nLoci]; if (arrayRawHo == NULL) memoryError(); for (countL = 0; countL < nLoci; countL++) { arrayRawHo[countL] = new float* [(nYear + 1)]; if (arrayRawHo == NULL) memoryError(); for (countY = 0; countY < (nYear + 1); countY++) { arrayRawHo[countL][countY] = new float [nIteration]; if (arrayRawHo == NULL) memoryError(); } } // He [nLoci][nYear][nIteration] float ***arrayRawHe; arrayRawHe = new float** [nLoci]; if (arrayRawHe == NULL) memoryError(); for (countL = 0; countL < nLoci; countL++) { arrayRawHe[countL] = new float* [(nYear + 1)]; if (arrayRawHe == NULL) memoryError(); for (countY = 0; countY < (nYear + 1); countY++) { arrayRawHe[countL][countY] = new float [nIteration]; if (arrayRawHe == NULL) memoryError(); } } // F [nLoci][nYear][nIteration] float ***arrayRawF; arrayRawF = new float** [nLoci]; if (arrayRawF == NULL) memoryError(); for (countL = 0; countL < nLoci; countL++) { arrayRawF[countL] = new float* [(nYear + 1)]; if (arrayRawF == NULL) memoryError(); for (countY = 0; countY < (nYear + 1); countY++) { arrayRawF[countL][countY] = new float [nIteration]; if (arrayRawF == NULL) memoryError(); } } // vector for list // allele pool for generating the pre-bottleneck population vector vectAllelePool; // population list // list of all individuals in the population // for selecting surviving individuals vector vectPopList; // start simulation // start of iteration loop for (countI = 0; countI < nIteration; countI++) { // generate the pre bottleneck population for (countL = 0; countL < nLoci; countL++) { // reset allele pool vector vectAllelePool.erase(vectAllelePool.begin(), vectAllelePool.end()); // fill allele pool vector with the number of allele counts from input data for (countA = 0; countA < maxAllele; countA++) { for (countAId = 0; countAId < arrayAlleleCount[countL][countA]; countAId++) { vectAllelePool.push_back(countA); } } // fill pre bottleneck population genotype array for (countN = 0; countN < nPopSize[0]; countN++) { for (countA = 0; countA < allelePerLocus; countA++) { // random shuffle allele pool random_shuffle(vectAllelePool.begin(), vectAllelePool.end()); // get allele from allele pool arrayGenotype[countN][countL][countA] = vectAllelePool[0]; } } } // end of countL loop for generating pre bottleneck population // stat of pre bottleneck population // call statDiploidMultiLocus function statDiploidMultiLocus(arrayGenotype, allelePerLocus, nPopSize[0], nLoci, arrayOA, arrayEA, arrayHo, arrayHe, arrayF); // write to raw data array for (countL = 0; countL < nLoci; countL++) { arrayRawOA[countL][0][countI] = arrayOA[countL]; arrayRawEA[countL][0][countI] = arrayEA[countL]; arrayRawHo[countL][0][countI] = arrayHo[countL]; arrayRawHe[countL][0][countI] = arrayHe[countL]; arrayRawF[countL][0][countI] = arrayF[countL]; } // fill age array for (countN = 0; countN < nPopSize[0]; countN++) { // age // when degreeOverlap = 0, if-test will always return true -> all individual start at age 0 // when degreeOverlap = 100, if-test will always return fales -> all individual start with random age value if (degreeOverlap < getRandomNumber(1, 100)) { arrayAge[countN] = 0; } else { arrayAge[countN] = getRandomNumber(0, (longevity - 1)); } } // start of year loop // each year loop has 2 major steps // 1. from previous year end to current year begin // generate survived population from the previous year // 2. from current year begin to current year end // generate breeding population of the current year // check age/replace old individuals reached longevity limit // geneerate new individuals if population growth occurred // // individual death events (old genotypes eliminated or replaced by new ones) can happen at both steps // death in step 1 is due to population decline (random death) // death in step 2 is due to reaching longevity limit (check age) // individual birth events always happen at step 2 // possible birth event1: old individuals reached longevity limit // possible birth event2: population growth occurred for (countY = 0; countY < nYear; countY++) { // check if asexual/monoecious or dioecious // if asexual/monoecious --> treat the population as a whole // else (dioecious) --> treat female and male separately if (reproSys >= 1 && reproSys <= 4) // asexual/monoecious, treat the population aas a whole { // step1: from previous year end to curent year begin // generate survived population (arrayGenotypeTemp) // check if population decline/stasis/growth // if decline -> get survivor list, copy survivor genotype and age // if stasis/growth -> copy all genotype and age from previous year end // check if population decline if (nPopSize[(countY+1)] < nPopSize[countY]) { nSurvivedTotal = nPopSize[(countY + 1)]; // fill population list vector for selecting surviving individuals for (countN = 0; countN < nPopSize[countY]; countN++) vectPopList.push_back(countN); // random suffle to generate survivor list random_shuffle(vectPopList.begin(), vectPopList.end()); // copy survivor genotype from arrayGenotype (previous year end) to arrayGenotypeTemp (current year begin) // copy survivor age (and increase by 1) for (countN = 0; countN < nSurvivedTotal; countN++) { for (countL = 0; countL < nLoci; countL++) { arrayGenotypeTemp[countN][countL][0] = arrayGenotype[(vectPopList[countN])][countL][0]; arrayGenotypeTemp[countN][countL][1] = arrayGenotype[(vectPopList[countN])][countL][1]; } arrayAgeTemp[countN] = (arrayAge[(vectPopList[countN])] + 1); } // end of copy genotype/age for loop // erase population list vector for next year vectPopList.erase(vectPopList.begin(), vectPopList.end()); } // end of if population decline else // population stasis or growth { nSurvivedTotal = nPopSize[countY]; // copy all genotype and age from previous year for (countN = 0; countN < nSurvivedTotal; countN++) { for (countL = 0; countL < nLoci; countL++) { arrayGenotypeTemp[countN][countL][0] = arrayGenotype[countN][countL][0]; arrayGenotypeTemp[countN][countL][1] = arrayGenotype[countN][countL][1]; } arrayAgeTemp[countN] = (arrayAge[countN] + 1); } } // end of population stasis or growth else // step2: from current year begin to current year end // generate breeding population // check individual age // if age reached longevity limit -> generate new individuals // check if population growth -> generate new individuals // copy genotype from arrayGenotypeTemp (current year begin) to arrayGenotype (current year end) // copy age from arrayAgeTemp to arrayAge // if population decline /stasis // -> total nPopSize[(countY+1)] individuals, all are from the previous year // if population growth // -> total nPopSize[(countY+1)] individuals, nPopSize[countY] are from the previous year // -> individuals from (nPopSize[countY]+1) to nPopSize[(countY+1)] are newborns of current year // generate breeding population // check if reached reproductive maturation, copy genotype to arrayParent nAdultTotal = 0; // reset number of adults for (countN = 0; countN < nSurvivedTotal; countN++) { if (arrayAgeTemp[countN] >= ageMaturation) { for (countL = 0; countL < nLoci; countL++) { arrayParent[nAdultTotal][countL][0] = arrayGenotypeTemp[countN][countL][0]; arrayParent[nAdultTotal][countL][1] = arrayGenotypeTemp[countN][countL][1]; } nAdultTotal++; } } // check survivor age // replace with new genotype if reached limit // else copy from arrayGenotypeTemp to arrayGenotype for (countN = 0; countN < nSurvivedTotal; countN++) { // check if reached longevity limit, generate new individual if (arrayAgeTemp[countN] >= longevity) { // reset age arrayAge[countN] = 0; // generate new genotype // total nAdultTotal breeders getNewDiMLGenotype(arrayParent, arrayNewGenotype, reproSys, nAdultTotal, nAdultFemale, parent1, parent2, nLoci); // copy genotype for (countL = 0; countL < nLoci; countL++) { arrayGenotype[countN][countL][0] = arrayNewGenotype[countL][0]; arrayGenotype[countN][countL][1] = arrayNewGenotype[countL][1]; } } // end of check age if else { // copy age arrayAge[countN] = arrayAgeTemp[countN]; // copy genotype (from current year begin to current year end) for (countL = 0; countL < nLoci; countL++) { arrayGenotype[countN][countL][0] = arrayGenotypeTemp[countN][countL][0]; arrayGenotype[countN][countL][1] = arrayGenotypeTemp[countN][countL][1]; } } } // end of check age for loop // if population growth occurs // individual from (nPopSize[countY]) to (nPopSize[(countY+1)] - 1) are newborns of current year for (countN = nPopSize[countY]; countN < nPopSize[(countY+1)]; countN++) { // reset age (newborn age = 0) arrayAge[countN] = 0; // generate new genotype // total nAdultTotal breeders getNewDiMLGenotype(arrayGenotypeTemp, arrayNewGenotype, reproSys, nAdultTotal, nAdultFemale, parent1, parent2, nLoci); // copy genotype for (countL = 0; countL < nLoci; countL++) { arrayGenotype[countN][countL][0] = arrayNewGenotype[countL][0]; arrayGenotype[countN][countL][1] = arrayNewGenotype[countL][1]; } } // end of population growth for loop } // end of asexual/monoecious if else // dioecious reproSys, treat female/male separately { // step1: from previous year end to curent year begin // generate breeding population (arrayGenotypeTemp) // check if nFemale/nMale decline/stasis/growth // if decline -> get survivor list, copy survivor genotype and age // if stasis/growth -> copy all genotype and age from previous year end // check if nFemale decline if (nFemale[(countY+1)] < nFemale[countY]) { // number of survived female = nFemale[(countY + 1)] nSurvivedFemale = nFemale[(countY + 1)]; // fill nFemale list vector for selecting surviving individuals for (countN = 0; countN < nFemale[countY]; countN++) vectPopList.push_back(countN); // random suffle to generate survivor list random_shuffle(vectPopList.begin(), vectPopList.end()); // copy survivor genotype from arrayGenotype (previous year end) to arrayGenotypeTemp (current year begin) // copy survivor age (and increase by 1) for (countN = 0; countN < nSurvivedFemale; countN++) { for (countL = 0; countL < nLoci; countL++) { arrayGenotypeTemp[countN][countL][0] = arrayGenotype[(vectPopList[countN])][countL][0]; arrayGenotypeTemp[countN][countL][1] = arrayGenotype[(vectPopList[countN])][countL][1]; } arrayAgeTemp[countN] = (arrayAge[(vectPopList[countN])] + 1); } // end of copy genotype/age for loop // erase population list vector for next year vectPopList.erase(vectPopList.begin(), vectPopList.end()); } // end of if nFemale decline else // nFemale stasis or growth { // number of survived female = nFemale[countY] nSurvivedFemale = nFemale[countY]; // copy all genotype and age from previous year for (countN = 0; countN < nSurvivedFemale; countN++) { for (countL = 0; countL < nLoci; countL++) { arrayGenotypeTemp[countN][countL][0] = arrayGenotype[countN][countL][0]; arrayGenotypeTemp[countN][countL][1] = arrayGenotype[countN][countL][1]; } arrayAgeTemp[countN] = (arrayAge[countN] + 1); } } // end of nFemale stasis or growth else // check if nMale decline if (nMale[(countY+1)] < nMale[countY]) { // number of survived male = nMale[(countY + 1)] nSurvivedMale = nMale[(countY + 1)]; // fill nMale list vector for selecting surviving individuals for (countN = 0; countN < nMale[countY]; countN++) vectPopList.push_back((countN + nFemale[countY])); // random suffle to generate survivor list random_shuffle(vectPopList.begin(), vectPopList.end()); // copy survivor genotype from arrayGenotype (previous year end) to arrayGenotypeTemp (current year begin) // copy survivor age (and increase by 1) for (countN = 0; countN < nSurvivedMale; countN++) { for (countL = 0; countL < nLoci; countL++) { arrayGenotypeTemp[(countN + nSurvivedFemale)][countL][0] = arrayGenotype[(vectPopList[countN])][countL][0]; arrayGenotypeTemp[(countN + nSurvivedFemale)][countL][1] = arrayGenotype[(vectPopList[countN])][countL][1]; } arrayAgeTemp[(countN + nSurvivedFemale)] = (arrayAge[(vectPopList[countN])] + 1); } // end of copy genotype/age for loop // erase population list vector for next year vectPopList.erase(vectPopList.begin(), vectPopList.end()); } // end of if nMale decline else // nMale stasis or growth { // number of survived male = nMale[countY] nSurvivedMale = nMale[countY]; // copy all genotype and age from previous year for (countN = 0; countN < nSurvivedMale; countN++) { for (countL = 0; countL < nLoci; countL++) { arrayGenotypeTemp[(countN + nSurvivedFemale)][countL][0] = arrayGenotype[(countN + nFemale[countY])][countL][0]; arrayGenotypeTemp[(countN + nSurvivedFemale)][countL][1] = arrayGenotype[(countN + nFemale[countY])][countL][1]; } arrayAgeTemp[(countN + nSurvivedFemale)] = (arrayAge[(countN + nFemale[countY])] + 1); } } // end of nMale stasis or growth else // step2: from current year begin to current year end // check individual age // if age reached longevity limit -> generate new individuals // check if population growth -> generate new individuals // copy genotype from arrayGenotypeTemp (current year begin) to arrayGenotype (current year end) // copy age from arrayAgeTemp to arrayAge // generate female/male newborns if needed // generate breeding population nAdultFemale = 0; for (countN = 0; countN < nSurvivedFemale; countN++) { if (arrayAgeTemp[countN] >= ageMaturation) { for (countL = 0; countL < nLoci; countL++) { arrayParent[nAdultFemale][countL][0] = arrayGenotypeTemp[countN][countL][0]; arrayParent[nAdultFemale][countL][1] = arrayGenotypeTemp[countN][countL][1]; } nAdultFemale++; } } nAdultMale = 0; for (countN = 0; countN < nSurvivedMale; countN++) { if (arrayAgeTemp[(countN + nSurvivedFemale)] >= ageMaturation) { for (countL = 0; countL < nLoci; countL++) { arrayParent[(nAdultMale + nAdultFemale)][countL][0] = arrayGenotypeTemp[(countN + nSurvivedFemale)][countL][0]; arrayParent[(nAdultMale + nAdultFemale)][countL][1] = arrayGenotypeTemp[(countN + nSurvivedFemale)][countL][1]; } nAdultMale++; } } nAdultTotal = (nAdultFemale + nAdultMale); // get parent id for diecious model with single reproducing male/pair // single reproducing male // get male id if (reproSys == 6) { parent2 = getRandomNumber(nAdultFemale, (nAdultTotal - 1)); } // single reproducing pair // get both parent ids if (reproSys == 7) { parent1 = getRandomNumber(0, (nAdultFemale - 1)); parent2 = getRandomNumber(nAdultFemale, (nAdultTotal - 1)); } // check survived female age // replace with new genotype if reached limit // else copy from arrayGenotypeTemp to arrayGenotype for (countN = 0; countN < nSurvivedFemale; countN++) { // check if reached longevity limit, generate new individual if (arrayAgeTemp[countN] >= longevity) { // reset age arrayAge[countN] = 0; // generate new genotype // total nAdultTotal breeders getNewDiMLGenotype(arrayGenotypeTemp, arrayNewGenotype, reproSys, nAdultTotal, nAdultFemale, parent1, parent2, nLoci); // copy genotype for (countL = 0; countL < nLoci; countL++) { arrayGenotype[countN][countL][0] = arrayNewGenotype[countL][0]; arrayGenotype[countN][countL][1] = arrayNewGenotype[countL][1]; } } // end of check age if else { // copy age arrayAge[countN] = arrayAgeTemp[countN]; // copy genotype (from current year begin to current year end) for (countL = 0; countL < nLoci; countL++) { arrayGenotype[countN][countL][0] = arrayGenotypeTemp[countN][countL][0]; arrayGenotype[countN][countL][1] = arrayGenotypeTemp[countN][countL][1]; } } } // individual from (nFemale[countY]) to (nFemale[(countY+1)] - 1) are newborns of current year for (countN = nFemale[countY]; countN < nFemale[(countY+1)]; countN++) { // reset age (newborn age = 0) arrayAge[countN] = 0; // generate new genotype // total nAdultTotal breeders getNewDiMLGenotype(arrayGenotypeTemp, arrayNewGenotype, reproSys, nAdultTotal, nAdultFemale, parent1, parent2, nLoci); // copy genotype for (countL = 0; countL < nLoci; countL++) { arrayGenotype[countN][countL][0] = arrayNewGenotype[countL][0]; arrayGenotype[countN][countL][1] = arrayNewGenotype[countL][1]; } } // end of nFemale growth for loop // check survived male age // replace with new genotype if reached limit // else copy from arrayGenotypeTemp to arrayGenotype for (countN = 0; countN < nSurvivedMale; countN++) { // check if reached longevity limit, generate new individual if (arrayAgeTemp[(countN + nSurvivedFemale)] >= longevity) { // reset age arrayAge[(countN + nFemale[(countY + 1)])] = 0; // generate new genotype // total nAdultTotal breeders getNewDiMLGenotype(arrayGenotypeTemp, arrayNewGenotype, reproSys, nAdultTotal, nAdultFemale, parent1, parent2, nLoci); // copy genotype for (countL = 0; countL < nLoci; countL++) { arrayGenotype[(countN + nFemale[(countY + 1)])][countL][0] = arrayNewGenotype[countL][0]; arrayGenotype[(countN + nFemale[(countY + 1)])][countL][1] = arrayNewGenotype[countL][1]; } } // end of check age if else { // copy age arrayAge[(countN + nFemale[(countY + 1)])] = arrayAgeTemp[(countN + nSurvivedFemale)]; // copy genotype (from current year begin to current year end) for (countL = 0; countL < nLoci; countL++) { arrayGenotype[(countN + nFemale[(countY + 1)])][countL][0] = arrayGenotypeTemp[(countN + nSurvivedFemale)][countL][0]; arrayGenotype[(countN + nFemale[(countY + 1)])][countL][1] = arrayGenotypeTemp[(countN + nSurvivedFemale)][countL][1]; } } } // end of check age for loop // if nMale growth occurs for (countN = 0; countN < (nMale[(countY+1)] - nMale[countY]); countN++) { // reset age (newborn age = 0) arrayAge[(countN + nFemale[(countY + 1)])] = 0; // generate new genotype // total nAdultTotal breeders getNewDiMLGenotype(arrayGenotypeTemp, arrayNewGenotype, reproSys, nAdultTotal, nAdultFemale, parent1, parent2, nLoci); // copy genotype for (countL = 0; countL < nLoci; countL++) { arrayGenotype[(countN + nFemale[(countY + 1)] + nSurvivedMale)][countL][0] = arrayNewGenotype[countL][0]; arrayGenotype[(countN + nFemale[(countY + 1)] + nSurvivedMale)][countL][1] = arrayNewGenotype[countL][1]; } } // end of nMale growth for loop } // end of dioecious reproSys else // stat of bottleneck population // call statDiploidMultiLocus function statDiploidMultiLocus(arrayGenotype, allelePerLocus, nPopSize[(countY+1)], nLoci, arrayOA, arrayEA, arrayHo, arrayHe, arrayF); // write to raw data array for (countL = 0; countL < nLoci; countL++) { arrayRawOA[countL][(countY + 1)][countI] = arrayOA[countL]; arrayRawEA[countL][(countY + 1)][countI] = arrayEA[countL]; arrayRawHo[countL][(countY + 1)][countI] = arrayHo[countL]; arrayRawHe[countL][(countY + 1)][countI] = arrayHe[countL]; arrayRawF[countL][(countY + 1)][countI] = arrayF[countL]; } } // end of year loop // check boolOutputRawGenotype if (boolOutputRawGenotype == 1) { // call writeDiSLRawGenotype function writeDiMLRawGenotype(rawGenotypeOutputFile, arrayGenotype, allelePerLocus, nPopSize[nYear], nLoci); } cout << '.'; // use dot to show iteration progress } // end of iteration loop // display message cout << endl; cout << "Simulation completed." << endl; cout << endl; cout << "Preparing output file..." << endl; cout << endl; // write simulation settings to output file // call writeVPSSimSetting function writeVPSSimSetting(outputFile, degreeOverlap, reproSys, longevity, ageMaturation, nYear, nIteration); // write summary output file writeDiMLSummary(outputFile, nLoci, arrayInputLocus, nYear, nIteration, arrayRawOA, arrayRawEA, arrayRawHo, arrayRawHe, arrayRawF); // call closeFile function // close input file closeInputFile(genotypeInputFile); closeInputFile(popSizeInput); // close output file closeOutputFile(outputFile); // check boolOutputRawGenotype if (boolOutputRawGenotype == 1) { closeOutputFile(rawGenotypeOutputFile); } } // end of simDiploidMultiLocusVPS function // end of source code file