set up two patches landscape, progress on #44

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<int> 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);
 	}
 }
 

NeutralStatsManager.cpp

@@ -39,7 +39,7 @@ NeutralStatsManager::NeutralStatsManager(const int& nbSampledPatches, const int
 void NeutralStatsManager::updateAllNeutralTables(Species* pSpecies, Landscape* pLandscape, set<int> 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);

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);

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<Patch*> patches(nbPatches);
+		set<int> 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<int> stgToSample = { 1 };
+
+		// Create genetic structure
+		const int genomeSz = 4;
+		Species* pSpecies = new Species();
+		pSpecies->setDemogr(createDefltDiploidDemogrParams());
+		pSpecies->setGeneticParameters(
+			set<int>{genomeSz - 1}, // single chromosome
+			genomeSz,
+			0.0, // no recombination
+			patchList,
+			indSampling,
+			stgToSample,
+			1
+		);
+		const set<int> 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<NeutralStatsManager>(nbPatches, nbLoci);
+		pNeutralStatistics->updateAllNeutralTables(
+			pSpecies, 
+			pLandscape, 
+			patchList
+		);
+		const int maxNbNeutralAlleles = static_cast<int>(maxAlleleVal);
+		pNeutralStatistics->calculateFstatWC(
+			patchList,
+			nbIndsPerPop,
+			nbLoci,
+			maxNbNeutralAlleles,
+			pSpecies,
+			pLandscape
+		);
+		assert(pNeutralStatistics->getFstWC() == 1.0);
 	}
 }