diff --git a/Community.cpp b/Community.cpp index b31a235..5f7a181 100644 --- a/Community.cpp +++ b/Community.cpp @@ -1845,7 +1845,7 @@ void Community::writePairwiseFSTFile(Species* pSpecies, const int yr, const int void Community::outNeutralGenetics(Species* pSpecies, int rep, int yr, int gen, bool fstat, bool perLocus, bool pairwise) { - const int nAlleles = (int)(pSpecies->getSpTrait(NEUTRAL)->getMutationParameters().find(MAX)->second); + const int maxNbNeutralAlleles = (int)(pSpecies->getSpTrait(NEUTRAL)->getMutationParameters().find(MAX)->second); const int nLoci = (int)pSpecies->getNPositionsForTrait(NEUTRAL); const set patchList = pSpecies->getSamplePatches(); int nInds = 0; @@ -1869,19 +1869,19 @@ void Community::outNeutralGenetics(Species* pSpecies, int rep, int yr, int gen, if (fstat) { pNeutralStatistics->calculateHo(patchList, nInds, nLoci, pSpecies, pLandscape); pNeutralStatistics->calcAllelicDiversityMetrics(patchList, nInds, pSpecies, pLandscape); - pNeutralStatistics->calculateFstatWC(patchList, nInds, nLoci, nAlleles, pSpecies, pLandscape); + pNeutralStatistics->calculateFstatWC(patchList, nInds, nLoci, maxNbNeutralAlleles, pSpecies, pLandscape); writeWCFstatFile(rep, yr, gen); } if (perLocus) { - pNeutralStatistics->calcPerLocusMeanSquaresFst(patchList, nInds, nLoci, nAlleles, pSpecies, pLandscape); + pNeutralStatistics->calcPerLocusMeanSquaresFst(patchList, nInds, nLoci, maxNbNeutralAlleles, pSpecies, pLandscape); pNeutralStatistics->calculatePerLocusHo(patchList, nInds, nLoci, pSpecies, pLandscape); - writeWCPerLocusFstatFile(pSpecies, yr, gen, nAlleles, nLoci, patchList); + writeWCPerLocusFstatFile(pSpecies, yr, gen, maxNbNeutralAlleles, nLoci, patchList); } if (pairwise) { pNeutralStatistics->calcPairwiseWeightedFst(patchList, nInds, nLoci, pSpecies, pLandscape); - writePairwiseFSTFile(pSpecies, yr, gen, nAlleles, nLoci, patchList); + writePairwiseFSTFile(pSpecies, yr, gen, maxNbNeutralAlleles, nLoci, patchList); } } diff --git a/NeutralStatsManager.cpp b/NeutralStatsManager.cpp index 307bce3..f64f023 100644 --- a/NeutralStatsManager.cpp +++ b/NeutralStatsManager.cpp @@ -39,7 +39,7 @@ NeutralStatsManager::NeutralStatsManager(const int& nbSampledPatches, const int void NeutralStatsManager::updateAllNeutralTables(Species* pSpecies, Landscape* pLandscape, set const& patchList) { const int nLoci = pSpecies->getNPositionsForTrait(NEUTRAL); - const int nAlleles = (int)pSpecies->getSpTrait(NEUTRAL)->getMutationParameters().find(MAX)->second; + const int maxNbNeutralAlleles = (int)pSpecies->getSpTrait(NEUTRAL)->getMutationParameters().find(MAX)->second; const int ploidy = pSpecies->isDiploid() ? 2 : 1; // Create / Update community-level NEUTRAL counts table @@ -48,7 +48,7 @@ void NeutralStatsManager::updateAllNeutralTables(Species* pSpecies, Landscape* p } else { // populate the tables with default values for (int thisLocus = 0; thisLocus < nLoci; thisLocus++) { - NeutralCountsTable newNeutralTbl = NeutralCountsTable(nAlleles); + NeutralCountsTable newNeutralTbl = NeutralCountsTable(maxNbNeutralAlleles); commNeutralCountTables.push_back(newNeutralTbl); } } @@ -67,7 +67,7 @@ void NeutralStatsManager::updateAllNeutralTables(Species* pSpecies, Landscape* p } // Add population-level counts to community-level counts for (int thisLocus = 0; thisLocus < nLoci; thisLocus++) { - for (int allele = 0; allele < nAlleles; allele++) { + for (int allele = 0; allele < maxNbNeutralAlleles; allele++) { if (pPop != 0) { patchAlleleCount = pPop->getAlleleTally(thisLocus, allele); diff --git a/NeutralStatsManager.h b/NeutralStatsManager.h index ee0fdf4..10e1594 100644 --- a/NeutralStatsManager.h +++ b/NeutralStatsManager.h @@ -50,7 +50,7 @@ struct PatchMatrix struct NeutralCountsTable { public: - NeutralCountsTable(int nAllele) : alleleTallies(nAllele), alleleFrequencies(nAllele), alleleHeterozygoteTallies(nAllele) {}; + NeutralCountsTable(int maxNbNeutralAlleles) : alleleTallies(maxNbNeutralAlleles), alleleFrequencies(maxNbNeutralAlleles), alleleHeterozygoteTallies(maxNbNeutralAlleles) {}; void reset() { fill(alleleTallies.begin(), alleleTallies.end(), 0); fill(alleleFrequencies.begin(), alleleFrequencies.end(), 0); diff --git a/unit_tests/testNeutralStats.cpp b/unit_tests/testNeutralStats.cpp index 9e1b5c5..604d504 100644 --- a/unit_tests/testNeutralStats.cpp +++ b/unit_tests/testNeutralStats.cpp @@ -174,7 +174,64 @@ void testNeutralStats() { // - per-locus Fst // - pairwise per-patch Fst { + // Create two-patches landscape + Landscape* pLandscape = new Landscape; + const int nbPatches = 2; + vector patches(nbPatches); + set patchList; + for (int i = 0; i < nbPatches; i++) { + patches[i] = pLandscape->newPatch(i); + Cell* pCell = new Cell(0, 0, (intptr)patches[i], 0); + patches[i]->addCell(pCell, 0, 0); + patchList.insert(patches[i]->getPatchNum()); + } + const string indSampling = "all"; + const set stgToSample = { 1 }; + + // Create genetic structure + const int genomeSz = 4; + Species* pSpecies = new Species(); + pSpecies->setDemogr(createDefltDiploidDemogrParams()); + pSpecies->setGeneticParameters( + set{genomeSz - 1}, // single chromosome + genomeSz, + 0.0, // no recombination + patchList, + indSampling, + stgToSample, + 1 + ); + const set genePositions = { 0, 1, 3 }; // arbitrary + const int nbLoci = genePositions.size(); + const bool isDiploid{ true }; + const float maxAlleleVal = 255; // highly unlikely that same value sampled twice + SpeciesTrait* spTr = createTestNeutralSpTrait(maxAlleleVal, genePositions, isDiploid); + pSpecies->addTrait(TraitType::NEUTRAL, *spTr); + // Initialise populations + const int nbIndsPerPop = 2; + for (auto p : patches) { + Population* pPop = new Population(pSpecies, p, nbIndsPerPop, 1); + pPop->sampleIndsWithoutReplacement(indSampling, stgToSample); + } + + // Calculate heterozygosity + auto pNeutralStatistics = make_unique(nbPatches, nbLoci); + pNeutralStatistics->updateAllNeutralTables( + pSpecies, + pLandscape, + patchList + ); + const int maxNbNeutralAlleles = static_cast(maxAlleleVal); + pNeutralStatistics->calculateFstatWC( + patchList, + nbIndsPerPop, + nbLoci, + maxNbNeutralAlleles, + pSpecies, + pLandscape + ); + assert(pNeutralStatistics->getFstWC() == 1.0); } }