SubCommunity.cpp

/*----------------------------------------------------------------------------
 *
 *	Copyright (C) 2026 Greta Bocedi, Stephen C.F. Palmer, Justin M.J. Travis, Anne-Kathleen Malchow, Roslyn Henry, Théo Pannetier, Jette Wolff, Damaris Zurell
 *
 *	This file is part of RangeShifter.
 *
 *	RangeShifter is free software: you can redistribute it and/or modify
 *	it under the terms of the GNU General Public License as published by
 *	the Free Software Foundation, either version 3 of the License, or
 *	(at your option) any later version.
 *
 *	RangeShifter is distributed in the hope that it will be useful,
 *	but WITHOUT ANY WARRANTY; without even the implied warranty of
 *	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 *	GNU General Public License for more details.
 *
 *	You should have received a copy of the GNU General Public License
 *	along with RangeShifter. If not, see <https://www.gnu.org/licenses/>.
 *
 --------------------------------------------------------------------------*/


 //---------------------------------------------------------------------------

#include "SubCommunity.h"
//---------------------------------------------------------------------------

ofstream outtraits;

//---------------------------------------------------------------------------

SubCommunity::SubCommunity(Patch* pPch, int num) {
	subCommNum = num;
	pPatch = pPch;
	// record the new sub-community no. in the patch
	pPatch->setSubComm(this);
	initial = false;
	occupancy = 0;
}

SubCommunity::~SubCommunity() {
	pPatch->setSubComm(nullptr);
	int npops = (int)popns.size();
	for (int i = 0; i < npops; i++) { // all populations
		delete popns[i];
	}
	popns.clear();
	if (occupancy != 0) delete[] occupancy;
}

int SubCommunity::getNum(void) { return subCommNum; }

Patch* SubCommunity::getPatch(void) { return pPatch; }

locn SubCommunity::getLocn(void) {
	locn loc = pPatch->getCellLocn(0);
	return loc;
}

void SubCommunity::setInitial(bool b) { initial = b; }

void SubCommunity::initialise(Landscape* pLandscape, Species* pSpecies)
{
	int ncells;
	landParams ppLand = pLandscape->getLandParams();
	initParams init = paramsInit->getInit();

	// determine size of initial population
	int nInds = 0;
	if (subCommNum == 0 // matrix patch
		|| !initial)   		// not in initial region or distribution
		nInds = 0;
	else {
		float k = pPatch->getK();
		if (k > 0.0) { // patch is currently suitable for this species
			switch (init.initDens) {
			case 0: // at carrying capacity
				nInds = (int)k;
				break;
			case 1: // at half carrying capacity
				nInds = (int)(k / 2.0);
				break;
			case 2: // specified no. per cell or density
				ncells = pPatch->getNCells();
				if (ppLand.patchModel) {
					nInds = (int)(init.indsHa * (float)(ncells * ppLand.resol * ppLand.resol) / 10000.0);
				}
				else {
					nInds = init.indsCell * ncells;
				}
				break;
			}
		}
		else nInds = 0;
	}

	// create new population only if it is non-zero or the matrix popn
	if (subCommNum == 0 || nInds > 0) {
		newPopn(pLandscape, pSpecies, pPatch, nInds);
	}

}

// initialise a specified individual
void SubCommunity::initialInd(Landscape* pLandscape, Species* pSpecies,
	Patch* pPatch, Cell* pCell, int ix)
{
	demogrParams dem = pSpecies->getDemogrParams();
	stageParams sstruct = pSpecies->getStageParams();
	emigRules emig = pSpecies->getEmigRules();
	transferRules trfr = pSpecies->getTransferRules();
	settleType sett = pSpecies->getSettle();
	short stg, age, repInt;
	Individual* pInd;
	float probmale;

	// create new population if not already in existence
	int npopns = (int)popns.size();
	if (npopns < 1) {
		newPopn(pLandscape, pSpecies, pPatch, 0);
	}

	// create new individual
	initInd iind = paramsInit->getInitInd(ix);
	if (dem.stageStruct) {
		stg = iind.stage; age = iind.age;   repInt = sstruct.repInterval;
	}
	else {
		age = stg = 1;   repInt = 0;
	}
	if (dem.repType == 0) {
		probmale = 0.0;
	}
	else {
		if (iind.sex == 1) probmale = 1.0; else probmale = 0.0;
	}
	pInd = new Individual(pSpecies, pCell, pPatch, stg, age, repInt, probmale, trfr.usesMovtProc, trfr.moveType);

	if (pSpecies->getNTraits() > 0) {
		// individual variation - set up genetics
		landData land = pLandscape->getLandData();
		pInd->setUpGenes(pSpecies, land.resol);
	}

	popns[0]->recruit(pInd);
}

// Create a new population, and return its address
Population* SubCommunity::newPopn(Landscape* pLandscape, Species* pSpecies,
	Patch* pPatch, int nInds)
{
	landParams land = pLandscape->getLandParams();
	int npopns = (int)popns.size();
	popns.push_back(new Population(pSpecies, pPatch, nInds, land.resol));
	return popns[npopns];
}

popStats SubCommunity::getPopStats(void) {
	popStats p, pop;
	p.pSpecies = 0; p.spNum = 0; p.nInds = p.nAdults = p.nNonJuvs = 0; p.breeding = false;
	p.pPatch = pPatch;
	// FOR SINGLE SPECIES IMPLEMENTATION, THERE IS ONLY ONE POPULATION IN THE PATCH
	int npops = (int)popns.size();
	for (int i = 0; i < npops; i++) { // all populations
		pop = popns[i]->getStats(pPatch->getDemoScaling());
		p.pSpecies = pop.pSpecies;
		p.spNum = pop.spNum;
		p.nInds += pop.nInds;
		p.nNonJuvs += pop.nNonJuvs;
		p.nAdults += pop.nAdults;
		p.breeding = pop.breeding;
	}
	return p;
}

void SubCommunity::resetPopns(void) {
	int npops = (int)popns.size();
	for (int i = 0; i < npops; i++) { // all populations
		delete popns[i];
	}
	popns.clear();
	// clear the list of populations in the corresponding patch
	pPatch->resetPopn();
}

void SubCommunity::resetPossSettlers(void) {
	if (subCommNum == 0) return; // not applicable in the matrix
	pPatch->resetPossSettlers();
}

// Extirpate all populations according to
// option 0 - random local extinction probability
// option 1 - local extinction probability gradient
// NB only applied for cell-based model
void SubCommunity::localExtinction(int option) {
	double pExtinct = 0.0;
	if (option == 0) {
		envStochParams env = paramsStoch->getStoch();
		if (env.localExt) pExtinct = env.locExtProb;
	}
	else {
		envGradParams grad = paramsGrad->getGradient();
		Cell* pCell = pPatch->getRandomCell(); // get only cell in the patch
		// extinction prob is complement of cell gradient value plus any non-zero prob at the optimum
		pExtinct = 1.0 - pCell->getEnvVal() + grad.extProbOpt;
		if (pExtinct > 1.0) pExtinct = 1.0;
	}
	if (pRandom->Bernoulli(pExtinct)) {
		int npops = (int)popns.size();
		for (int i = 0; i < npops; i++) { // all populations
			popns[i]->extirpate();
		}
	}
}

// Action in event of patch becoming unsuitable owing to landscape change
void SubCommunity::patchChange(void) {
	if (subCommNum == 0) return; // no reproduction in the matrix
	Species* pSpecies;
	float localK = 0.0;
	int npops = (int)popns.size();
	// THE FOLLOWING MAY BE MORE EFFICIENT WHILST THERE IS ONLY ONE SPECIES ...
	if (npops < 1) return;
	localK = pPatch->getK();
	if (localK <= 0.0) { // patch in dynamic landscape has become unsuitable
		for (int i = 0; i < npops; i++) { // all populations
			pSpecies = popns[i]->getSpecies();
			demogrParams dem = pSpecies->getDemogrParams();
			if (dem.stageStruct) {
				stageParams sstruct = pSpecies->getStageParams();
				if (sstruct.disperseOnLoss) popns[i]->allEmigrate();
				else popns[i]->extirpate();
			}
			else { // non-stage-structured species is destroyed
				popns[i]->extirpate();
			}
		}
	}
}

void SubCommunity::reproduction(int resol, float epsGlobal, short rasterType, bool patchModel)
{
	if (subCommNum == 0) return; // no reproduction in the matrix
	float localK, envval;
	std::vector <float> localDemoScaling;
	Cell* pCell;
	envGradParams grad = paramsGrad->getGradient();
	envStochParams env = paramsStoch->getStoch();

	int npops = (int)popns.size();
	// THE FOLLOWING MAY BE MORE EFFICIENT WHILST THERE IS ONLY ONE SPECIES ...
	if (npops < 1) return;

	localK = pPatch->getK();
	localDemoScaling = pPatch->getDemoScaling();
	if (localK > 0.0) {
		if (patchModel) {
			envval = 1.0; // environmental gradient is currently not applied for patch-based model
		}
		else { // cell-based model
			if (grad.gradient && grad.gradType == 2)
			{ // gradient in fecundity
				Cell* pCell = pPatch->getRandomCell(); // locate the only cell in the patch
				envval = pCell->getEnvVal();
			}
			else envval = 1.0;
		}
		if (env.stoch && !env.inK) { // stochasticity in fecundity
			if (env.local) {
				if (!patchModel) { // only permitted for cell-based model
					pCell = pPatch->getRandomCell();
					if (pCell != 0) envval += pCell->getEps();
				}
			}
			else { // global stochasticity
				envval += epsGlobal;
			}
		}
		for (int i = 0; i < npops; i++) { // all populations
			popns[i]->reproduction(localK, envval, resol, localDemoScaling);
			popns[i]->fledge();
		}
	}
}

void SubCommunity::emigration(void)
{
	if (subCommNum == 0) return; // no emigration from the matrix
	int npops = static_cast<int>(popns.size());
	if (npops < 1) return;
	float localK = pPatch->getK();
	// NOTE that even if K is zero, it could have been >0 in previous time-step, and there
	// might be emigrants if there is non-juvenile emigration
	for (int i = 0; i < npops; i++) { // all populations
		popns[i]->emigration(localK);
	}
}

// Remove emigrants from their natal patch and add to a map of vectors
void SubCommunity::recruitDispersers(std::map<Species*,std::vector<Individual *>>& disperserPool) {
	if (subCommNum == 0) return; // no dispersal initiation in the matrix
	popStats pop;
	disperser disp;

	int npops = (int)popns.size();
	for (int i = 0; i < npops; i++) { // all populations
		pop = popns[i]->getStats(pPatch->getDemoScaling());
		Species* pSpecies = popns[i]->getSpecies();
		for (int j = 0; j < pop.nInds; j++) {
			disp = popns[i]->extractDisperser(j);
			if (disp.yes) { // disperser - has already been removed from natal population
				disperserPool[pSpecies].push_back(disp.pInd);
			}
		}
		// remove pointers to emigrants
		popns[i]->clean();
	}
}

// Add all individuals in the matrix to the disperser pool
void SubCommunity::disperseMatrix(std::map<Species*,std::vector<Individual *>> &inds_map) {
	if (subCommNum != 0) return;
	popStats pop;

	int npops = (int)popns.size();
	for (int i = 0; i < npops; i++) {
		pop = popns[i]->getStats(pPatch->getDemoScaling());
		Species* pSpecies = popns[i]->getSpecies();
#pragma omp for schedule(static)
		for (int j = 0; j < pop.nInds; j++) {
			Individual *pInd = popns[i]->extractIndividual(j);
			inds_map[pSpecies].push_back(pInd);
}
#pragma omp single
		popns[i]->clean();
	}

}

// Add an individual into the local population of its species in the patch
void SubCommunity::recruit(Individual* pInd, Species* pSpecies) {
	int npops = (int)popns.size();
	for (int i = 0; i < npops; i++) { // all populations
		if (pSpecies == popns[i]->getSpecies()) {
			popns[i]->recruit(pInd);
		}
	}
}

// Add individuals into the local population of their species in the patch
void SubCommunity::recruitMany(std::vector<Individual*>& inds, Species* pSpecies) {
	int npops = (int)popns.size();
	for (int i = 0; i < npops; i++) { // all populations
		if (pSpecies == popns[i]->getSpecies()) {
			popns[i]->recruitMany(inds);
		}
	}
}

// Transfer through the matrix - run for a per-species map of vectors of individuals
int SubCommunity::resolveTransfer(std::map<Species*,vector<Individual*>>& dispersingInds, Landscape* pLandscape, short landIx)
{
	int nbStillDispersing = 0;
	int disperser;
	Patch* pPatch = nullptr;
	Cell* pCell = nullptr;
	simParams sim = paramsSim->getSim();

	for (auto & it : dispersingInds) { // all species

		Species* const& pSpecies = it.first;
		short reptype = pSpecies->getRepType();
		transferRules trfr = pSpecies->getTransferRules();

		vector<Individual*>& inds = it.second;
		
		// each individual takes one step
		// for dispersal by kernel, this should be the only step taken
		for (auto& pInd : inds) {
			if (trfr.usesMovtProc) {
				disperser = pInd->moveStep(pLandscape, pSpecies, landIx, sim.absorbing);
			}
			else {
				disperser = pInd->moveKernel(pLandscape, pSpecies, sim.absorbing);
			}
			nbStillDispersing += disperser;
			if (disperser) {
				if (reptype > 0)
				{ // sexual species - record as potential settler in new patch
					if (pInd->getStatus() == 2)
					{ // disperser has found a patch
						pCell = pInd->getCurrCell();
						pPatch = pCell->getPatch();
						if (pPatch != nullptr) { // not no-data area
							pPatch->incrPossSettler(pSpecies, pInd->getSex());
						}
					}
				}
			}
		}
	}
	return nbStillDispersing;
}


// Determine whether there is a potential mate present in a patch which a potential
// settler has reached
bool SubCommunity::matePresent(Species* pSpecies, Cell* pCell, short othersex)
{
	Patch* pPatch;
	Population* pNewPopn;
	int popsize = 0;
	bool matefound = false;

	pPatch = pCell->getPatch();
	if (pPatch != nullptr) {
		if (pPatch->getPatchNum() > 0) { // not the matrix patch
			if (pPatch->getK() > 0.0)
			{ // suitable
				pNewPopn = pPatch->getPopn(pSpecies);
				if (pNewPopn != nullptr) {
					const stageParams sstruct = pSpecies->getStageParams();
					// count members of other sex already resident in the patch
					for (int stg = 0; stg < sstruct.nStages; stg++) {
						popsize += pNewPopn->getNbInds(stg, othersex);
					}
				}
				if (popsize < 1) {
					// add any potential settlers of the other sex
					popsize += pPatch->getPossSettlers(pSpecies, othersex);
				}
			}
		}
	}
	if (popsize > 0) matefound = true;
	return matefound;
}

// Transfer is run for populations in the matrix only
#if RS_RCPP // included also SEASONAL
int SubCommunity::resolveSettlement(std::map<Species*, vector<Individual*>>& dispersingInds, Landscape* pLandscape, short nextseason)
#else
int SubCommunity::resolveSettlement(std::map<Species*, vector<Individual*>>& dispersingInds, Landscape* pLandscape)
#endif
{
	int nbStillDispersing = 0;
	short othersex;
	bool mateOK, densdepOK;
	int patchnum;
	double localK, popsize, settprob;
	Patch* pPatch = 0;
	Cell* pCell = 0;
	indStats ind;
	Population* pNewPopn = 0;
	locn newloc, nbrloc;

	landData ppLand = pLandscape->getLandData();
	settleRules sett;
	settleTraits settDD;
	settlePatch settle;
	simParams sim = paramsSim->getSim();

	for (auto& it : dispersingInds) { // all species
		Species* const& pSpecies = it.first;
		transferRules trfr = pSpecies->getTransferRules();
		settleType settletype = pSpecies->getSettle();

		vector<Individual*>& inds = it.second;
		
		// each individual which has reached a potential patch decides whether to settle
		for (auto& pInd : inds) {
			ind = pInd->getStats();
			if (ind.sex == 0) othersex = 1; else othersex = 0;
			if (settletype.stgDep) {
				if (settletype.sexDep) sett = pSpecies->getSettRules(ind.stage, ind.sex);
				else sett = pSpecies->getSettRules(ind.stage, 0);
			}
			else {
				if (settletype.sexDep) sett = pSpecies->getSettRules(0, ind.sex);
				else sett = pSpecies->getSettRules(0, 0);
			}
			if (ind.status == 2)
			{ // awaiting settlement
				pCell = pInd->getCurrCell();
				if (pCell == 0) {
					// this condition can occur in a patch-based model at the time of a dynamic landscape
					// change when there is a range restriction in place, since a patch can straddle the
					// range restriction and an individual forced to disperse upon patch removal could
					// start its trajectory beyond the boundary of the restrictyed range - such a model is 
					// not good practice, but the condition must be handled by killing the individual conceerned
					ind.status = 6;
				}
				else {
					mateOK = false;
					if (sett.findMate) {
						// determine whether at least one individual of the opposite sex is present in the
						// new population
						if (matePresent(pSpecies, pCell, othersex)) mateOK = true;
					}
					else { // no requirement to find a mate
						mateOK = true;
					}

					densdepOK = false;
					settle = pInd->getSettPatch();
					if (sett.densDep)
					{
						pPatch = pCell->getPatch();
						if (pPatch != nullptr) { // not no-data area
							if (settle.settleStatus == 0
								|| settle.pSettPatch != pPatch)
								// note: second condition allows for having moved from one patch to another
								// adjacent one
							{
								// determine whether settlement occurs in the (new) patch
								localK = (double)pPatch->getK();
								pNewPopn = pPatch->getPopn(pSpecies);
								if (pNewPopn == nullptr) { // population has not been set up in the new patch
									popsize = 0.0;
								}
								else {
									popsize = (double)pNewPopn->getNbInds();
								}
								if (localK > 0.0) {
									// make settlement decision
									if (settletype.indVar) settDD = pInd->getIndSettTraits();
#if RS_RCPP
									else settDD = pSpecies->getSpSettTraits(ind.stage, ind.sex);
#else
									else {
										if (settletype.sexDep) {
											if (settletype.stgDep)
												settDD = pSpecies->getSpSettTraits(ind.stage, ind.sex);
											else
												settDD = pSpecies->getSpSettTraits(0, ind.sex);
										}
										else {
											if (settletype.stgDep)
												settDD = pSpecies->getSpSettTraits(ind.stage, 0);
											else
												settDD = pSpecies->getSpSettTraits(0, 0);
										}
									}
#endif //RS_RCPP
									settprob = settDD.s0 /
										(1.0 + exp(-(popsize / localK - (double)settDD.beta) * (double)settDD.alpha));

									if (pRandom->Bernoulli(settprob)) { // settlement allowed
										densdepOK = true;
										settle.settleStatus = 2;
									}
									else { // settlement procluded
										settle.settleStatus = 1;
									}
									settle.pSettPatch = pPatch;
								}
								pInd->setSettPatch(settle);
							}
							else {
								if (settle.settleStatus == 2) { // previously allowed to settle
									densdepOK = true;
								}
							}
						}
					}
					else { // no density-dependent settlement
						densdepOK = true;
						settle.settleStatus = 2;
						settle.pSettPatch = pPatch;
						pInd->setSettPatch(settle);
					}

					if (mateOK && densdepOK) { // can recruit to patch
						ind.status = 4;
						nbStillDispersing--;
					}
					else { // does not recruit
						if (trfr.usesMovtProc) {
							ind.status = 1; // continue dispersing, unless ...
							// ... maximum steps has been exceeded
							pathSteps steps = pInd->getSteps();
							settleSteps settsteps = pSpecies->getSteps(ind.stage, ind.sex);
							if (steps.year >= settsteps.maxStepsYr) {
								ind.status = 3; // waits until next year
							}
							if (steps.total >= settsteps.maxSteps) {
								ind.status = 6; // dies
							}
						}
						else { // dispersal kernel
							if (sett.wait) {
								ind.status = 3; // wait until next dispersal event
							}
							else {
								ind.status = 6; // (dies unless a neighbouring cell is suitable)
							}
							nbStillDispersing--;
						}
					}
				}

				pInd->setStatus(ind.status);
			}
#if RS_RCPP
			// write each individuals current movement step and status to paths file
			if (trfr.usesMovtProc && sim.outPaths) {
				if (nextseason >= sim.outStartPaths && nextseason % sim.outIntPaths == 0) {
					pInd->outMovePath(nextseason);
				}
			}
#endif

			if (!trfr.usesMovtProc && sett.go2nbrLocn && (ind.status == 3 || ind.status == 6))
			{
				// for kernel-based transfer only ...
				// determine whether recruitment to a neighbouring cell is possible

				pCell = pInd->getCurrCell();
				newloc = pCell->getLocn();
				vector <Cell*> nbrlist;
				for (int dx = -1; dx < 2; dx++) {
					for (int dy = -1; dy < 2; dy++) {
						if (dx != 0 || dy != 0) { //cell is not the current cell
							nbrloc.x = newloc.x + dx; nbrloc.y = newloc.y + dy;
							if (nbrloc.x >= 0 && nbrloc.x <= ppLand.maxX
								&& nbrloc.y >= 0 && nbrloc.y <= ppLand.maxY) { // within landscape
								// add to list of potential neighbouring cells if suitable, etc.
								pCell = pLandscape->findCell(nbrloc.x, nbrloc.y);
								if (pCell != 0) { // not no-data area
									pPatch = pCell->getPatch();
									if (pPatch != nullptr) { // not no-data area
										patchnum = pPatch->getPatchNum();
										if (patchnum > 0 && pPatch != pInd->getNatalPatch())
										{ // not the matrix or natal patch
											if (pPatch->getK() > 0.0)
											{ // suitable
												if (sett.findMate) {
													if (matePresent(pSpecies, pCell, othersex)) nbrlist.push_back(pCell);
												}
												else
													nbrlist.push_back(pCell);
											}
										}
									}
								}
							}
						}
					}
				}
				int listsize = (int)nbrlist.size();
				if (listsize > 0) { // there is at least one suitable neighbouring cell
					if (listsize == 1) {
						pInd->moveto(nbrlist[0]);
					}
					else { // select at random from the list
						int rrr = pRandom->IRandom(0, listsize - 1);
						pInd->moveto(nbrlist[rrr]);
					}
				}
				// else list empty - do nothing - individual retains its current location and status
			}

		} // loop through inds

	} // loop through disperser map

	return nbStillDispersing;
}

//---------------------------------------------------------------------------

// Remove emigrants from the vectors map and transfer to sub-community
// in which their destination co-ordinates fall

void SubCommunity::completeDispersal(std::map<Species*,vector<Individual*>>& inds_map, Landscape* pLandscape, bool connect)
{
	Population* pPop;
	Patch* pPrevPatch;
	Patch* pNewPatch;
	Cell* pPrevCell;
	SubCommunity* pSubComm;

	for (auto & it : inds_map) { // all species
		Species* const& pSpecies = it.first;
		vector<Individual*>& inds = it.second;
		for (Individual*& pInd : inds) {
			indStats ind = pInd->getStats();
			bool settled = ind.status == 4 || ind.status == 5;
			if (settled) {
				// find new patch
				pNewPatch = pInd->getCurrCell()->getPatch();
				// find population within the patch (if there is one)
				{
#ifdef _OPENMP
					const std::unique_lock<std::mutex> lock = pNewPatch->lockPopns();
#endif // _OPENMP
					pPop = pNewPatch->getPopn(pSpecies);
					if (pPop == 0) { // settler is the first in a previously uninhabited patch
						// create a new population in the corresponding sub-community
						pSubComm = pNewPatch->getSubComm();
						pPop = pSubComm->newPopn(pLandscape, pSpecies, pNewPatch, 0);
					}
				}
				pPop->recruit(pInd);
				if (connect) { // increment connectivity totals
					int newpatch = pNewPatch->getSeqNum();
					pPrevCell = pInd->getPrevCell();
					pPrevPatch = pPrevCell->getPatch();
					if (pPrevPatch != nullptr) {
						int prevpatch = pPrevPatch->getSeqNum();
						pLandscape->incrConnectMatrix(prevpatch, newpatch);
					}
				}
				pInd = nullptr;
			}
			else { // for group dispersal only
			}
		}
		// remove settled individuals
		inds.erase(std::remove(inds.begin(), inds.end(), (Individual *)nullptr), inds.end());
	}
}

//---------------------------------------------------------------------------

void SubCommunity::survival0(short option0, short option1)
{
	int npops = (int)popns.size();
    std::vector <float> localDemoScaling;
		localDemoScaling = pPatch->getDemoScaling();
	float localK = pPatch->getK();
	for (int i = 0; i < npops; i++) { // all populations
			popns[i]->survival0(localK, option0, option1, localDemoScaling);
	}
}

void SubCommunity::survival1()
{
	int npops = (int)popns.size();
	for (int i = 0; i < npops; i++) { // all populations
		popns[i]->survival1();
	}
}

void SubCommunity::survival(short part, short option0, short option1
) {
	if (part == 0) {
		return survival0(option0, option1);
}
	else {
		return survival1();
	}
}

void SubCommunity::ageIncrement(void) {
	int npops = (int)popns.size();
	for (int i = 0; i < npops; i++) { // all populations
		popns[i]->ageIncrement();
	}
}

// Find the population of a given species in a given patch
Population* SubCommunity::findPop(Species* pSp, Patch* pPch) {

	Population* pPop = 0;
	popStats pop;
	int npops = (int)popns.size();

	for (int i = 0; i < npops; i++) { // all populations
		pop = popns[i]->getStats(pPatch->getDemoScaling());
		if (pop.pSpecies == pSp && pop.pPatch == pPch) { // population located
			pPop = popns[i];
			break;
		}
		else pPop = 0;
	}
	return pPop;
}

//---------------------------------------------------------------------------

void SubCommunity::createOccupancy(int nrows) {
	if (occupancy != 0) deleteOccupancy();
	if (nrows > 0) {
		occupancy = new int[nrows];
		for (int i = 0; i < nrows; i++) occupancy[i] = 0;
	}
}

void SubCommunity::updateOccupancy(int row) {
	popStats pop;
	int npops = (int)popns.size();
	for (int i = 0; i < npops; i++) {
		pop = popns[i]->getStats(pPatch->getDemoScaling());
		if (pop.nInds > 0 && pop.breeding) {
			occupancy[row]++;
			i = npops;
		}
	}
}

int SubCommunity::getOccupancy(int row) {
	if (row >= 0) return occupancy[row];
	else return 0;
}

void SubCommunity::deleteOccupancy(void) {
	delete[] occupancy;
	occupancy = 0;
}

//---------------------------------------------------------------------------
// Close population file
bool SubCommunity::outPopFinishLandscape()
{
	bool fileOK;
	Population* pPop;

	// as all populations may have been deleted, set up a dummy one
	// species is not necessary
	pPop = new Population();
	fileOK = pPop->outPopFinishLandscape();
	delete pPop;
	return fileOK;
}

//---------------------------------------------------------------------------
// Open population file and write header record
bool SubCommunity::outPopStartLandscape(Landscape* pLandscape, Species* pSpecies)
{
	bool fileOK;
	Population* pPop;
	landParams land = pLandscape->getLandParams();

	// as no population has yet been created, set up a dummy one
	// species is necessary, as columns depend on stage and sex structure
	pPop = new Population(pSpecies, pPatch, 0, land.resol);
	fileOK = pPop->outPopStartLandscape(land.landNum, land.patchModel);
	delete pPop;
	return fileOK;
}

// Write records to population file
void SubCommunity::outPop(Landscape* pLandscape, int rep, int yr, int gen)
{
	landParams land = pLandscape->getLandParams();
	envGradParams grad = paramsGrad->getGradient();
	envStochParams env = paramsStoch->getStoch();
	bool writeEnv = false;
	bool gradK = false;
	if (grad.gradient) {
		writeEnv = true;
		if (grad.gradType == 1) gradK = true; // ... carrying capacity
	}
	if (env.stoch) writeEnv = true;

	// generate output for each population within the sub-community (patch)
	// provided that the patch is suitable (i.e. non-zero carrying capacity)
	// or the population is above zero (possible if there is stochasticity or a moving gradient)
	// or it is the matrix patch in a patch-based model
	int npops = (int)popns.size();
	int patchnum;
	Cell* pCell;
	float localK;
	float eps = 0.0;
	if (env.stoch) {
		if (env.local) {
			pCell = pPatch->getRandomCell();
			if (pCell != 0) eps = pCell->getEps();
		}
		else {
			eps = pLandscape->getGlobalStoch(yr);
		}
	}

	patchnum = pPatch->getPatchNum();
	for (int i = 0; i < npops; i++) { // all populations
		localK = pPatch->getK();
		if (localK > 0.0 || (land.patchModel && patchnum == 0)) {
			popns[i]->outPopulation(rep, yr, gen, eps, land.patchModel, writeEnv, gradK);
		}
		else {
			if (popns[i]->getNbInds() > 0) {
				popns[i]->outPopulation(rep, yr, gen, eps, land.patchModel, writeEnv, gradK);
			}
		}
	}
}

// Close individuals file
void SubCommunity::outIndsFinishReplicate() {
	// as all populations have been deleted, set up a dummy one
	Population* pPop = new Population();
	pPop->outIndsFinishReplicate();
	delete pPop;
		return;
}

// Open individuals file and write header record
void SubCommunity::outIndsStartReplicate(Landscape* pLandscape, int rep, int landNr) {
	landParams ppLand = pLandscape->getLandParams();
	popns[0]->outIndsStartReplicate(rep, landNr, ppLand.patchModel);
}

// Write records to individuals file
void SubCommunity::outIndividuals(Landscape* pLandscape, int rep, int yr, int gen) {
	// generate output for each population within the sub-community (patch)
	int npops = (int)popns.size();
	for (int i = 0; i < npops; i++) { // all populations
		popns[i]->outIndividual(pLandscape, rep, yr, gen, pPatch->getPatchNum());
	}
}


// Population size of a specified stage
int SubCommunity::getNbInds(int stage) const {
	int popsize = 0;
	int npops = (int)popns.size();
	for (int i = 0; i < npops; i++) { // all populations
		popsize += popns[i]->getNbInds(stage);
	}
	return popsize;
}

// Close traits file
bool SubCommunity::outTraitsFinishLandscape()
{
	if (outtraits.is_open()) outtraits.close();
	outtraits.clear();
	return true;
}

// Open traits file and write header record
bool SubCommunity::outTraitsStartLandscape(Landscape* pLandscape, Species* pSpecies, int landNr)
{
	landParams land = pLandscape->getLandParams();
	string name;
	emigRules emig = pSpecies->getEmigRules();
	transferRules trfr = pSpecies->getTransferRules();
	settleType sett = pSpecies->getSettle();
	simParams sim = paramsSim->getSim();

	string DirOut = paramsSim->getDir(2);
	if (sim.batchMode) {
		if (land.patchModel) {
			name = DirOut
				+ "Batch" + to_string(sim.batchNum) + "_"
				+ "Sim" + to_string(sim.simulation) + "_Land" + to_string(landNr)
				+ "_TraitsXpatch.txt";
		}
		else {
			name = DirOut
				+ "Batch" + to_string(sim.batchNum) + "_"
				+ "Sim" + to_string(sim.simulation) + "_Land" + to_string(landNr)
				+ "_TraitsXcell.txt";
		}
	}
	else {
		if (land.patchModel) {
			name = DirOut + "Sim" + to_string(sim.simulation) + "_TraitsXpatch.txt";
		}
		else {
			name = DirOut + "Sim" + to_string(sim.simulation) + "_TraitsXcell.txt";
		}
	}
	outtraits.open(name.c_str());

	outtraits << "Rep\tYear\tRepSeason";
	if (land.patchModel) outtraits << "\tPatchID";
	else
		outtraits << "\tx\ty";

	if (emig.indVar) {
		if (emig.sexDep) {
			if (emig.densDep) {
				outtraits << "\tF_meanD0\tF_stdD0\tM_meanD0\tM_stdD0";
				outtraits << "\tF_meanAlpha\tF_stdAlpha\tM_meanAlpha\tM_stdAlpha";
				outtraits << "\tF_meanBeta\tF_stdBeta\tM_meanBeta\tM_stdBeta";
			}
			else {
				outtraits << "\tF_meanEP\tF_stdEP\tM_meanEP\tM_stdEP";
			}
		}
		else {
			if (emig.densDep) {
				outtraits << "\tmeanD0\tstdD0\tmeanAlpha\tstdAlpha";
				outtraits << "\tmeanBeta\tstdBeta";
			}
			else {
				outtraits << "\tmeanEP\tstdEP";
			}
		}
	}
	if (trfr.indVar) {
		if (trfr.usesMovtProc) {
			if (trfr.moveType == 1) {
				outtraits << "\tmeanDP\tstdDP\tmeanGB\tstdGB";
				outtraits << "\tmeanAlphaDB\tstdAlphaDB\tmeanBetaDB\tstdBetaDB";
			}
			if (trfr.moveType == 2) {
				outtraits << "\tmeanStepLength\tstdStepLength\tmeanRho\tstdRho";
			}
		}
		else {
			if (trfr.sexDep) {
				outtraits << "\tF_mean_distI\tF_std_distI\tM_mean_distI\tM_std_distI";
				if (trfr.twinKern)
					outtraits << "\tF_mean_distII\tF_std_distII\tM_mean_distII\tM_std_distII"
					<< "\tF_meanPfirstKernel\tF_stdPfirstKernel"
					<< "\tM_meanPfirstKernel\tM_stdPfirstKernel";
			}
			else {
				outtraits << "\tmean_distI\tstd_distI";
				if (trfr.twinKern)
					outtraits << "\tmean_distII\tstd_distII\tmeanPfirstKernel\tstdPfirstKernel";
			}
		}
	}
	if (sett.indVar) {
		if (sett.sexDep) {
			outtraits << "\tF_meanS0\tF_stdS0\tM_meanS0\tM_stdS0";
			outtraits << "\tF_meanAlphaS\tF_stdAlphaS\tM_meanAlphaS\tM_stdAlphaS";
			outtraits << "\tF_meanBetaS\tF_stdBetaS\tM_meanBetaS\tM_stdBetaS";
		}
		else {
			outtraits << "\tmeanS0\tstdS0";
			outtraits << "\tmeanAlphaS\tstdAlphaS";
			outtraits << "\tmeanBetaS\tstdBetaS";
		}
	}
	if (pSpecies->getNbGenLoadTraits() > 0) {
		if (pSpecies->getDemogrParams().repType > 0) {
			outtraits << "\tF_meanGenFitness\tF_stdGenFitness\tM_meanGenFitness\tM_stdGenFitness";
		}
		else {
			outtraits << "\tmeanGenFitness\tstdGenFitness";
		}
	}

	outtraits << endl;

	return outtraits.is_open();
}

// Write records to traits file and return aggregated sums
traitsums SubCommunity::outTraits(Landscape* pLandscape, int rep, int yr, int gen, bool commlevel)
{
	landParams land = pLandscape->getLandParams();
	simParams sim = paramsSim->getSim();
	traitsums ts, indTraitsSums;
	for (int i = 0; i < gMaxNbSexes; i++) {
		ts.ninds[i] = 0;
		ts.sumD0[i] = ts.ssqD0[i] = 0.0;
		ts.sumAlpha[i] = ts.ssqAlpha[i] = 0.0; ts.sumBeta[i] = ts.ssqBeta[i] = 0.0;
		ts.sumDist1[i] = ts.ssqDist1[i] = 0.0; ts.sumDist2[i] = ts.ssqDist2[i] = 0.0;
		ts.sumProp1[i] = ts.ssqProp1[i] = 0.0;
		ts.sumDP[i] = ts.ssqDP[i] = 0.0;
		ts.sumGB[i] = ts.ssqGB[i] = 0.0;
		ts.sumAlphaDB[i] = ts.ssqAlphaDB[i] = 0.0;
		ts.sumBetaDB[i] = ts.ssqBetaDB[i] = 0.0;
		ts.sumStepL[i] = ts.ssqStepL[i] = 0.0; ts.sumRho[i] = ts.ssqRho[i] = 0.0;
		ts.sumS0[i] = ts.ssqS0[i] = 0.0;
		ts.sumAlphaS[i] = ts.ssqAlphaS[i] = 0.0; ts.sumBetaS[i] = ts.ssqBetaS[i] = 0.0;
		ts.sumGeneticFitness[i] = ts.ssqGeneticFitness[i] = 0.0;
	}

	// generate output for each population within the sub-community (patch)
	// provided that the patch is suitable (i.e. non-zero carrying capacity)
	int npops = (int)popns.size();
	Species* pSpecies;

	for (int iPop = 0; iPop < npops; iPop++) { // all populations

		if (pPatch->getK() > 0.0 && popns[iPop]->getNbInds() > 0) {
			pSpecies = popns[iPop]->getSpecies();
			demogrParams dem = pSpecies->getDemogrParams();
			emigRules emig = pSpecies->getEmigRules();
			transferRules trfr = pSpecies->getTransferRules();
			settleType sett = pSpecies->getSettle();
			indTraitsSums = popns[iPop]->getIndTraitsSums(pSpecies);

			// Sum over populations
			for (int iSex = 0; iSex < gMaxNbSexes; iSex++) {
				ts.ninds[iSex] += indTraitsSums.ninds[iSex];
				ts.sumD0[iSex] += indTraitsSums.sumD0[iSex];
				ts.ssqD0[iSex] += indTraitsSums.ssqD0[iSex];
				ts.sumAlpha[iSex] += indTraitsSums.sumAlpha[iSex];
				ts.ssqAlpha[iSex] += indTraitsSums.ssqAlpha[iSex];
				ts.sumBeta[iSex] += indTraitsSums.sumBeta[iSex];
				ts.ssqBeta[iSex] += indTraitsSums.ssqBeta[iSex];
				ts.sumDist1[iSex] += indTraitsSums.sumDist1[iSex];
				ts.ssqDist1[iSex] += indTraitsSums.ssqDist1[iSex];
				ts.sumDist2[iSex] += indTraitsSums.sumDist2[iSex];
				ts.ssqDist2[iSex] += indTraitsSums.ssqDist2[iSex];
				ts.sumProp1[iSex] += indTraitsSums.sumProp1[iSex];
				ts.ssqProp1[iSex] += indTraitsSums.ssqProp1[iSex];
				ts.sumDP[iSex] += indTraitsSums.sumDP[iSex];
				ts.ssqDP[iSex] += indTraitsSums.ssqDP[iSex];
				ts.sumGB[iSex] += indTraitsSums.sumGB[iSex];
				ts.ssqGB[iSex] += indTraitsSums.ssqGB[iSex];
				ts.sumAlphaDB[iSex] += indTraitsSums.sumAlphaDB[iSex];
				ts.ssqAlphaDB[iSex] += indTraitsSums.ssqAlphaDB[iSex];
				ts.sumBetaDB[iSex] += indTraitsSums.sumBetaDB[iSex];
				ts.ssqBetaDB[iSex] += indTraitsSums.ssqBetaDB[iSex];
				ts.sumStepL[iSex] += indTraitsSums.sumStepL[iSex];
				ts.ssqStepL[iSex] += indTraitsSums.ssqStepL[iSex];
				ts.sumRho[iSex] += indTraitsSums.sumRho[iSex];
				ts.ssqRho[iSex] += indTraitsSums.ssqRho[iSex];
				ts.sumS0[iSex] += indTraitsSums.sumS0[iSex];
				ts.ssqS0[iSex] += indTraitsSums.ssqS0[iSex];
				ts.sumAlphaS[iSex] += indTraitsSums.sumAlphaS[iSex];
				ts.ssqAlphaS[iSex] += indTraitsSums.ssqAlphaS[iSex];
				ts.sumBetaS[iSex] += indTraitsSums.sumBetaS[iSex];
				ts.ssqBetaS[iSex] += indTraitsSums.ssqBetaS[iSex];
				ts.sumGeneticFitness[iSex] += indTraitsSums.sumGeneticFitness[iSex];
				ts.ssqGeneticFitness[iSex] += indTraitsSums.ssqGeneticFitness[iSex];
			}

			// Produce trait-per-cell output
			if (sim.outTraitsCells && yr % sim.outIntTraitCell == 0
				&& !commlevel) {

				outtraits << rep << "\t" << yr << "\t" << gen;
				if (land.patchModel) {
					outtraits << "\t" << pPatch->getPatchNum();
				}
				else {
					locn loc = pPatch->getCellLocn(0);
					outtraits << "\t" << loc.x << "\t" << loc.y;
				}

			if (emig.indVar) {
					if (emig.sexDep) {
						vector<double> mnD0(2, 0.0), mnAlpha(2, 0.0), mnBeta(2, 0.0), sdD0(2, 0.0), sdAlpha(2, 0.0), sdBeta(2, 0.0);
						for (int sex = 0; sex < gMaxNbSexes; sex++) {

							double popsize = static_cast<double>(indTraitsSums.ninds[sex]);

					if (popsize > 0) {

								mnD0[sex] = indTraitsSums.sumD0[sex] / popsize;
								mnAlpha[sex] = indTraitsSums.sumAlpha[sex] / popsize;
								mnBeta[sex] = indTraitsSums.sumBeta[sex] / popsize;

						if (popsize > 1) {
									sdD0[sex] = indTraitsSums.ssqD0[sex] / popsize - mnD0[sex] * mnD0[sex];
									if (sdD0[sex] > 0.0) sdD0[sex] = sqrt(sdD0[sex]); else sdD0[sex] = 0.0;
									sdAlpha[sex] = indTraitsSums.ssqAlpha[sex] / popsize - mnAlpha[sex] * mnAlpha[sex];
									if (sdAlpha[sex] > 0.0) sdAlpha[sex] = sqrt(sdAlpha[sex]); else sdAlpha[sex] = 0.0;
									sdBeta[sex] = indTraitsSums.ssqBeta[sex] / popsize - mnBeta[sex] * mnBeta[sex];
									if (sdBeta[sex] > 0.0) sdBeta[sex] = sqrt(sdBeta[sex]); else sdBeta[sex] = 0.0;
						}
						else {
									sdD0[sex] = sdAlpha[sex] = sdBeta[sex] = 0.0;
						}
					}
				}
						outtraits << "\t" << mnD0[0] << "\t" << sdD0[0];
						outtraits << "\t" << mnD0[1] << "\t" << sdD0[1];
						if (emig.densDep) {
							outtraits << "\t" << mnAlpha[0] << "\t" << sdAlpha[0];
							outtraits << "\t" << mnAlpha[1] << "\t" << sdAlpha[1];
							outtraits << "\t" << mnBeta[0] << "\t" << sdBeta[0];
							outtraits << "\t" << mnBeta[1] << "\t" << sdBeta[1];
						}
					}
					else { // not sex-dependent
						double mnD0 = 0, mnAlpha = 0, mnBeta = 0, popsize = 0;
						double sdD0 = 0, sdAlpha = 0, sdBeta = 0;
						for (int sex = 0; sex < gMaxNbSexes; sex++) {
							mnD0 += indTraitsSums.sumD0[sex];
							mnAlpha += indTraitsSums.sumAlpha[sex];
							mnBeta += indTraitsSums.sumBeta[sex];
							popsize += indTraitsSums.ninds[sex];
							sdD0 += indTraitsSums.ssqD0[sex];
							sdAlpha += indTraitsSums.ssqAlpha[sex];
							sdBeta += indTraitsSums.ssqBeta[sex];
						}
						mnD0 /= popsize;
						mnAlpha /= popsize;
						mnBeta /= popsize;
						if (popsize > 1) {
							sdD0 = sdD0 / popsize - mnD0 * mnD0;
							sdAlpha = sdAlpha / popsize - mnAlpha * mnAlpha;
							sdBeta = sdBeta / popsize - mnBeta * mnBeta;
							sdD0 = sdD0 == 0.0 ? 0.0 : sqrt(sdD0);
							sdAlpha = sdAlpha == 0.0 ? 0.0 : sqrt(sdAlpha);
							sdBeta = sdBeta == 0.0 ? 0.0 : sqrt(sdBeta);
						}
						else {
							sdD0 = 0.0;
							sdAlpha = 0.0;
							sdBeta = 0.0;
						}
						
						outtraits << "\t" << mnD0 << "\t" << sdD0;
						if (emig.densDep) {
							outtraits << "\t" << mnAlpha << "\t" << sdAlpha;
							outtraits << "\t" << mnBeta << "\t" << sdBeta;
						}
					}
				}

			if (trfr.indVar) {

					if (trfr.usesMovtProc) { // not sex-dependent
						if (trfr.moveType == 1) { // SMS
							double mnDP = 0.0, mnGB = 0.0, mnAlphaDB = 0.0, mnBetaDB = 0.0;
							double sdDP = 0.0, sdGB = 0.0, sdAlphaDB = 0.0, sdBetaDB = 0.0;
							double popsize = 0.0;
							for (int sex = 0; sex < gMaxNbSexes; sex++) {
								mnDP += indTraitsSums.sumDP[sex];
								mnGB += indTraitsSums.sumGB[sex];
								mnAlphaDB += indTraitsSums.sumAlphaDB[sex];
								mnBetaDB+= indTraitsSums.sumBetaDB[sex];
								popsize += indTraitsSums.ninds[sex];
								sdDP += indTraitsSums.ssqDP[sex];
								sdGB += indTraitsSums.ssqGB[sex];
								sdAlphaDB += indTraitsSums.ssqAlphaDB[sex];
								sdBetaDB += indTraitsSums.ssqBetaDB[sex];
				}
							mnDP /= popsize;
							mnGB /= popsize;
							mnAlphaDB /= popsize;
							mnBetaDB /= popsize;
							if (popsize > 1) {
								sdDP = sdDP / popsize - mnDP * mnDP;
								sdGB = sdGB / popsize - mnGB * mnGB;
								sdAlphaDB = sdAlphaDB / popsize - mnAlphaDB * mnAlphaDB;
								sdBetaDB = sdBetaDB / popsize - mnBetaDB * mnBetaDB;
								sdDP = sdDP == 0.0 ? 0.0 : sqrt(sdDP);
								sdGB = sdGB == 0.0 ? 0.0 : sqrt(sdGB);
								sdAlphaDB = sdAlphaDB == 0.0 ? 0.0 : sqrt(sdAlphaDB);
								sdBetaDB = sdBetaDB == 0.0 ? 0.0 : sqrt(sdBetaDB);
							}
				else {
								sdDP = 0.0;
								sdGB = 0.0;
								sdAlphaDB = 0.0;
								sdBetaDB = 0.0;
					}
							outtraits << "\t" << mnDP << "\t" << sdDP;
							outtraits << "\t" << mnGB << "\t" << sdGB;
							outtraits << "\t" << mnAlphaDB << "\t" << sdAlphaDB;
							outtraits << "\t" << mnBetaDB << "\t" << sdBetaDB;
						}
						if (trfr.moveType == 2) { // CRW
							double mnStepL = 0.0, mnRho = 0.0, sdStepL = 0.0, sdRho = 0.0;
							double popsize = 0.0;
							for (int sex = 0; sex < gMaxNbSexes; sex++) {
								mnStepL += indTraitsSums.sumStepL[sex];
								mnRho += indTraitsSums.sumRho[sex];
								popsize += indTraitsSums.ninds[sex];
								sdStepL += indTraitsSums.ssqStepL[sex];
								sdRho += indTraitsSums.ssqRho[sex];
							}
							mnStepL /= popsize;
							mnRho /= popsize;
							if (popsize > 1) {
								sdStepL = sdStepL / popsize - mnStepL * mnStepL;
								sdRho = sdRho / popsize - mnRho * mnRho;
								sdStepL = sdStepL == 0.0 ? 0.0 : sqrt(sdStepL);
								sdRho = sdRho == 0.0 ? 0.0 : sqrt(sdRho);
							}
					else {
								sdStepL = 0.0;
								sdRho = 0.0;
					}
							outtraits << "\t" << mnStepL << "\t" << sdStepL;
							outtraits << "\t" << mnRho << "\t" << sdRho;
				}
					}
					else { // kernels
						if (trfr.sexDep) {

							vector<double> mnDist1(2, 0.0), mnDist2(2, 0.0), mnProp1(2, 0.0), mnStepL(2, 0.0), mnRho(2, 0.0);
							vector<double> sdDist1(2, 0.0), sdDist2(2, 0.0), sdProp1(2, 0.0), sdStepL(2, 0.0), sdRho(2, 0.0);

							for (int sex = 0; sex < gMaxNbSexes; sex++) {

					// individuals may have been counted by sex if there was
					// sex dependency in another dispersal phase
								double popsize = static_cast<double>(indTraitsSums.ninds[sex]);

					if (popsize > 0) {
									mnDist1[sex] = indTraitsSums.sumDist1[sex] / popsize;
									mnDist2[sex] = indTraitsSums.sumDist2[sex] / popsize;
									mnProp1[sex] = indTraitsSums.sumProp1[sex] / popsize;
						if (popsize > 1) {
										sdDist1[sex] = indTraitsSums.ssqDist1[sex] / popsize - mnDist1[sex] * mnDist1[sex];
										if (sdDist1[sex] > 0.0) sdDist1[sex] = sqrt(sdDist1[sex]); else sdDist1[sex] = 0.0;
										sdDist2[sex] = indTraitsSums.ssqDist2[sex] / popsize - mnDist2[sex] * mnDist2[sex];
										if (sdDist2[sex] > 0.0) sdDist2[sex] = sqrt(sdDist2[sex]); else sdDist2[sex] = 0.0;
										sdProp1[sex] = indTraitsSums.ssqProp1[sex] / popsize - mnProp1[sex] * mnProp1[sex];
										if (sdProp1[sex] > 0.0) sdProp1[sex] = sqrt(sdProp1[sex]); else sdProp1[sex] = 0.0;
										sdStepL[sex] = indTraitsSums.ssqStepL[sex] / popsize - mnStepL[sex] * mnStepL[sex];
						}
					}
				}
							outtraits << "\t" << mnDist1[0] << "\t" << sdDist1[0];
							outtraits << "\t" << mnDist1[1] << "\t" << sdDist1[1];
							if (trfr.twinKern) {
								outtraits << "\t" << mnDist2[0] << "\t" << sdDist2[0];
								outtraits << "\t" << mnDist2[1] << "\t" << sdDist2[1];
								outtraits << "\t" << mnProp1[0] << "\t" << sdProp1[0];
								outtraits << "\t" << mnProp1[1] << "\t" << sdProp1[1];
							}
						}
						else { // sex-independent
							double mnDist1 = 0.0, mnDist2 = 0.0, mnProp1 = 0.0, popsize = 0.0;
							double sdDist1 = 0.0, sdDist2 = 0.0, sdProp1 = 0.0;
							for (int sex = 0; sex < gMaxNbSexes; sex++) {
								mnDist1 += indTraitsSums.sumDist1[sex];
								mnDist2 += indTraitsSums.sumDist2[sex];
								mnProp1 += indTraitsSums.sumProp1[sex];
								popsize += indTraitsSums.ninds[sex];
								sdDist1 += indTraitsSums.ssqDist1[sex];
								sdDist2 += indTraitsSums.ssqDist2[sex];
								sdProp1 += indTraitsSums.ssqProp1[sex];
							}
							mnDist1 /= popsize;
							mnDist2 /= popsize;
							mnProp1 /= popsize;
							if (popsize > 1) {
								sdDist1 = sdDist1 / popsize - mnDist1 * mnDist1;
								sdDist2 = sdDist2 / popsize - mnDist2 * mnDist2;
								sdProp1 = sdProp1 / popsize - mnProp1 * mnProp1;
								sdDist1 = sdDist1 == 0.0 ? 0.0 : sqrt(sdDist1);
								sdDist2 = sdDist2 == 0.0 ? 0.0 : sqrt(sdDist2);
								sdProp1 = sdProp1 == 0.0 ? 0.0 : sqrt(sdProp1);
						}
							else {
								sdDist1 = 0.0;
								sdDist2 = 0.0;
								sdProp1 = 0.0;
					}
							outtraits << "\t" << mnDist1 << "\t" << sdDist1;
							if (trfr.twinKern) {
								outtraits << "\t" << mnDist2 << "\t" << sdDist2;
								outtraits << "\t" << mnProp1 << "\t" << sdProp1;
				}
			}
					}
				}

			if (sett.indVar) {

					if (sett.sexDep) {

						vector<double> mnS0(2, 0.0), mnAlpha(2, 0.0), mnBeta(2, 0.0), sdS0(2, 0.0), sdAlpha(2, 0.0), sdBeta(2, 0.0);

						for (int sex = 0; sex < gMaxNbSexes; sex++) {

					// individuals may have been counted by sex if there was
					// sex dependency in another dispersal phase
							double popsize = static_cast<double>(indTraitsSums.ninds[sex]);

					if (popsize > 0) {

								mnS0[sex] = indTraitsSums.sumS0[sex] / popsize;
								mnAlpha[sex] = indTraitsSums.sumAlphaS[sex] / popsize;
								mnBeta[sex] = indTraitsSums.sumBetaS[sex] / popsize;

						if (popsize > 1) {
									sdS0[sex] = indTraitsSums.ssqS0[sex] / popsize - mnS0[sex] * mnS0[sex];
									if (sdS0[sex] > 0.0) sdS0[sex] = sqrt(sdS0[sex]); else sdS0[sex] = 0.0;
									sdAlpha[sex] = indTraitsSums.ssqAlphaS[sex] / popsize - mnAlpha[sex] * mnAlpha[sex];
									if (sdAlpha[sex] > 0.0) sdAlpha[sex] = sqrt(sdAlpha[sex]); else sdAlpha[sex] = 0.0;
									sdBeta[sex] = indTraitsSums.ssqBetaS[sex] / popsize - mnBeta[sex] * mnBeta[sex];
									if (sdBeta[sex] > 0.0) sdBeta[sex] = sqrt(sdBeta[sex]); else sdBeta[sex] = 0.0;
						}
						else {
									sdS0[sex] = sdAlpha[sex] = sdBeta[sex] = 0.0;
						}
					}
				}

						outtraits << "\t" << mnS0[0] << "\t" << sdS0[0];
						outtraits << "\t" << mnS0[1] << "\t" << sdS0[1];
						outtraits << "\t" << mnAlpha[0] << "\t" << sdAlpha[0];
						outtraits << "\t" << mnAlpha[1] << "\t" << sdAlpha[1];
						outtraits << "\t" << mnBeta[0] << "\t" << sdBeta[0];
						outtraits << "\t" << mnBeta[1] << "\t" << sdBeta[1];
					}
					else { // sex-independent
						double mnS0 = 0, mnAlpha = 0, mnBeta = 0, popsize = 0;
						double sdS0 = 0, sdAlpha = 0, sdBeta = 0;
						for (int sex = 0; sex < gMaxNbSexes; sex++) {
							mnS0 += indTraitsSums.sumS0[sex];
							mnAlpha += indTraitsSums.sumAlphaS[sex];
							mnBeta += indTraitsSums.sumBetaS[sex];
							popsize += indTraitsSums.ninds[sex];
							sdS0 += indTraitsSums.ssqS0[sex];
							sdAlpha += indTraitsSums.ssqAlphaS[sex];
							sdBeta += indTraitsSums.ssqBetaS[sex];
					}
						mnS0 /= popsize;
						mnAlpha /= popsize;
						mnBeta /= popsize;
						if (popsize > 1) {
							sdS0 = sdS0 / popsize - mnS0 * mnS0;
							sdAlpha = sdAlpha / popsize - mnAlpha * mnAlpha;
							sdBeta = sdBeta / popsize - mnBeta * mnBeta;
							sdS0 = sdS0 == 0.0 ? 0.0 : sqrt(sdS0);
							sdAlpha = sdAlpha == 0.0 ? 0.0 : sqrt(sdAlpha);
							sdBeta = sdBeta == 0.0 ? 0.0 : sqrt(sdBeta);
				}
						else {
							sdS0 = 0.0;
							sdAlpha = 0.0;
							sdBeta = 0.0;
			}
						outtraits << "\t" << mnS0 << "\t" << sdS0;
						outtraits << "\t" << mnAlpha << "\t" << sdAlpha;
						outtraits << "\t" << mnBeta << "\t" << sdBeta;
					}
				}

				// Genetic load
				if (pSpecies->getNbGenLoadTraits() > 0) {

					if (pSpecies->getDemogrParams().repType > 0) { // sexual model
						vector<double> mnGenFitness(2, 0.0), sdGenFitness(2, 0.0);

						for (int sex = 0; sex < gMaxNbSexes; sex++) {
							double popsize = static_cast<double>(indTraitsSums.ninds[sex]);
							mnGenFitness[sex] = indTraitsSums.sumGeneticFitness[sex]
								/ popsize;
							if (popsize > 1) {
								sdGenFitness[sex] = indTraitsSums.ssqGeneticFitness[sex]
									/ popsize - mnGenFitness[sex]
									* mnGenFitness[sex];
								if (sdGenFitness[sex] > 0.0)
									sdGenFitness[sex] = sqrt(sdGenFitness[sex]);
								else sdGenFitness[sex] = 0.0;
			}
							else {
								sdGenFitness[sex] = 0.0;
		}
	}
						outtraits << "\t" << mnGenFitness[0] << "\t" << sdGenFitness[0];
						outtraits << "\t" << mnGenFitness[1] << "\t" << sdGenFitness[1];
					}
					else { // asexual

						double mnGenFitness = 0.0, popsize = 0.0, sdGenFitness= 0.0;
						for (int sex = 0; sex < gMaxNbSexes; sex++) {
							mnGenFitness += indTraitsSums.sumGeneticFitness[sex];
							popsize += indTraitsSums.ninds[sex];
							sdGenFitness += indTraitsSums.ssqGeneticFitness[sex];
						}
						mnGenFitness /= popsize;
						if (popsize > 1) {
							sdGenFitness = sdGenFitness / popsize - mnGenFitness * mnGenFitness;
							sdGenFitness = sdGenFitness == 0.0 ? 0.0 : sqrt(sdGenFitness);
						}
						else sdGenFitness = 0.0;
						outtraits << "\t" << mnGenFitness << "\t" << sdGenFitness;
					}
				}
				outtraits << endl;
			} // end trait-per-cell output

		} 
	} // end population loop
	return ts;
}

//---------------------------------------------------------------------------

//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------