diff --git a/.gitignore b/.gitignore
index 9bd72d2..c4a4ac5 100644
--- a/.gitignore
+++ b/.gitignore
@@ -15,6 +15,7 @@ build-Release/
*.workspace
*.mk
*.tags
+/build/
# Hidden source
/RangeShiftR/src/.*
diff --git a/RandomCheck.h b/Allele.h
similarity index 58%
rename from RandomCheck.h
rename to Allele.h
index b4992d2..9acfc97 100644
--- a/RandomCheck.h
+++ b/Allele.h
@@ -1,40 +1,36 @@
+#ifndef ALLELEH
+#define ALLELEH
+
/*----------------------------------------------------------------------------
- *
- * Copyright (C) 2020 Greta Bocedi, Stephen C.F. Palmer, Justin M.J. Travis, Anne-Kathleen Malchow, Damaris Zurell
- *
+ *
+ * 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 .
- *
+ *
+ * File Created by Roslyn Henry March 2023. Code adapted from NEMO (https://nemo2.sourceforge.io/)
--------------------------------------------------------------------------*/
-
-
-//---------------------------------------------------------------------------
-
-#ifndef RandomCheckH
-#define RandomCheckH
-
-#include
-using namespace std;
-
-#include "Parameters.h"
-#include "RSrandom.h"
-
-void randomCheck(void);
-extern paramSim *paramsSim;
-extern RSrandom *pRandom;
-//---------------------------------------------------------------------------
+class Allele {
+ const float value;
+ const float dominance;
+public:
+ Allele(float alleleValue, float alleleDominance) : value(alleleValue), dominance(alleleDominance) { }
+ ~Allele() {}
+ float getAlleleValue() const { return value; };
+ float getDominanceCoef() const { return dominance; };
+};
#endif
diff --git a/CMakeLists.txt b/CMakeLists.txt
index c14a3e3..06869ac 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,30 +1,35 @@
# Config file for compilation with CMake
-if (NOT batchmode) # that is, RScore as a standalone
+if(NOT batchmode) # that is, RScore as a standalone
cmake_minimum_required(VERSION 3.10)
# set the project name and version
- project(RScore VERSION 2.1.0)
+ project(RScore VERSION 3.0.0)
# specify the C++ standard
- set(CMAKE_CXX_STANDARD 17)
+ set(CMAKE_CXX_STANDARD 20)
set(CMAKE_CXX_STANDARD_REQUIRED True)
- add_executable(RScore Main.cpp Species.cpp Cell.cpp Community.cpp FractalGenerator.cpp Genome.cpp Individual.cpp Landscape.cpp Model.cpp Parameters.cpp Patch.cpp Population.cpp RandomCheck.cpp RSrandom.cpp SubCommunity.cpp Utils.cpp)
+
+ add_executable(RScore Main.cpp Species.cpp Cell.cpp Community.cpp FractalGenerator.cpp GeneticFitnessTrait.cpp Individual.cpp Landscape.cpp Management.cpp Model.cpp NeutralStatsManager.cpp Parameters.cpp Patch.cpp Population.cpp DispersalTrait.cpp RSrandom.cpp NeutralTrait.cpp SpeciesTrait.cpp SubCommunity.cpp Utils.cpp "unit_tests/testIndividual.cpp" "unit_tests/testNeutralStats.cpp" "unit_tests/testPopulation.cpp")
+
+ # turn on unit tests
+ add_compile_definitions("UNIT_TESTS")
else() # that is, RScore compiled as library within RangeShifter_batch
- add_library(RScore Species.cpp Cell.cpp Community.cpp FractalGenerator.cpp Genome.cpp Individual.cpp Landscape.cpp Model.cpp Parameters.cpp Patch.cpp Population.cpp RandomCheck.cpp RSrandom.cpp SubCommunity.cpp Utils.cpp)
+
+ add_library(RScore Species.cpp Cell.cpp Community.cpp FractalGenerator.cpp GeneticFitnessTrait.cpp Individual.cpp Landscape.cpp Management.cpp Model.cpp NeutralStatsManager.cpp Parameters.cpp Patch.cpp Population.cpp DispersalTrait.cpp RSrandom.cpp NeutralTrait.cpp SpeciesTrait.cpp SubCommunity.cpp Utils.cpp)
endif()
-# pass config definitions to compiler
-target_compile_definitions(RScore PRIVATE RSWIN64)
+if(OMP)
+ find_package(OpenMP COMPONENTS CXX)
+ if(OpenMP_CXX_FOUND)
+ target_link_libraries(RScore PUBLIC OpenMP::OpenMP_CXX)
+ endif()
+endif()
# enable LINUX_CLUSTER macro on Linux + macOS
if(CMAKE_SYSTEM_NAME STREQUAL "Linux" OR CMAKE_SYSTEM_NAME STREQUAL "Darwin")
add_compile_definitions("LINUX_CLUSTER")
endif()
-# Debug Mode by default, unless "release" is passed
-if(NOT DEFINED release)
- add_compile_definitions(RSDEBUG)
-endif()
if(NOT batchmode)
target_include_directories(RScore PUBLIC "${PROJECT_BINARY_DIR}")
-endif()
\ No newline at end of file
+endif()
diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md
index e1b62ff..fa25488 100644
--- a/CONTRIBUTING.md
+++ b/CONTRIBUTING.md
@@ -17,7 +17,7 @@ RangeShifter is distributed with three user interfaces, each living in their own
All three share the same source code for the core simulation (i.e., the actual model), which lives in this repo (RScore). Each of the interfaces keeps a copy of this core code in a subfolder called RScore, kept in sync with the RScore repo via a git subtree (see Git subtree usage section).
-⚠️ If you wish to propose a change to one of the interfaces, please do so in the corresponding repo: [RangeShifter batch mode](https://github.com/RangeShifter/RangeShifter_batch_dev), [RangeShiftR package](https://github.com/RangeShifter/RangeShiftR-package-dev).
+⚠ If you wish to propose a change to one of the interfaces, please do so in the corresponding repo: [RangeShifter batch mode](https://github.com/RangeShifter/RangeShifter_batch_dev), [RangeShiftR package](https://github.com/RangeShifter/RangeShiftR-package-dev).
*The RangeShifter GUI is currently being rewritten, and is not open source yet.
@@ -26,7 +26,7 @@ All three share the same source code for the core simulation (i.e., the actual m
#### Maintainers
- [@JetteReeg](https://github.com/JetteReeg): RScore repo and lead in R package
-- [@TheoPannetier](https://github.com/TheoPannetier): RScore repo and lead in batch mode
+
Maintainers are responsible for coordinating development efforts and ensuring that RangeShifter keeps building continuously.
@@ -61,12 +61,17 @@ In the meantime, we encourage contributors to work in small and frequent commits
Any changes regarding the RangeShifter core code should be done in this repository and can afterwards be synced with all interfaces using the git subtree feature (see [Git subtree](https://github.com/RangeShifter/RScore/tree/main?tab=readme-ov-file#usage-git-subtree) section in the README).
+
#### Bugs
-To report a bug, please [open an issue](https://github.com/RangeShifter/RangeShiftR-package-dev/issues/new), using the Bug Report template.
-Please do check if a related issue has already open on one of the other interfaces ([here](https://github.com/RangeShifter/RangeShifter_batch-dev/issues) for the batch interface or [here](https://github.com/RangeShifter/RangeShiftR-package-dev) for the R package interface).
+To report a bug, please [open an issue](https://github.com/RangeShifter/RangeShiftR-package/issues/new), using the Bug Report template.
+Please do check if a related issue has already open on one of the other interfaces ([here](https://github.com/RangeShifter/RangeShifter_batch/issues) for the batch interface or [here](https://github.com/RangeShifter/RangeShiftR-package) for the R package interface).
+
To propose a bug fix (thank you!!), please create and work on your own branch or fork, from either `main` or `develop` (preferred), and open a pull request when your fix is ready to be merged into the original branch.
+As a prerequisite for merging, please ensure that your version passes status check (that is, RScore can still build, and all unit tests are still satisfied).
+This can be seen in the Actions panel for every commit and at the bottom of the pull request.
+
Maintainers will review the pull request, possibly request changes, and eventually integrate the bug fix into RScore, and update the subtrees to bring the fix to all interfaces.
#### New features
diff --git a/Cell.cpp b/Cell.cpp
index 1c0f937..73eda84 100644
--- a/Cell.cpp
+++ b/Cell.cpp
@@ -1,25 +1,25 @@
/*----------------------------------------------------------------------------
- *
- * Copyright (C) 2020 Greta Bocedi, Stephen C.F. Palmer, Justin M.J. Travis, Anne-Kathleen Malchow, Damaris Zurell
- *
+ *
+ * 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 .
- *
+ *
--------------------------------------------------------------------------*/
-
-
+
+
//---------------------------------------------------------------------------
#include "Cell.h"
@@ -30,101 +30,82 @@
// Cell functions
-Cell::Cell(int xx,int yy,intptr patch,int hab)
+Cell::Cell(int xx,int yy,Patch *patch,int hab)
{
-x = xx; y = yy;
-pPatch = patch;
-envVal = 1.0; // default - no effect of any gradient
-envDev = eps = 0.0;
-habIxx.push_back(hab);
-#if RSDEBUG
-//DebugGUI(("Cell::Cell(): this=" + Int2Str((int)this)
-// + " x=" + Int2Str(x) + " y=" + Int2Str(y)
-// + " habIndex=" + Int2Str(habIndex)
-//).c_str());
-#endif
-visits = 0;
-smsData = 0;
+ x = xx; y = yy;
+ pPatch = patch;
+ envVal = 1.0; // default - no effect of any gradient
+ envDev = eps = 0.0;
+ habIxx.push_back(hab);
+ visits = 0;
+ smsData = 0;
}
-Cell::Cell(int xx,int yy,intptr patch,float hab)
+Cell::Cell(int xx,int yy,Patch *patch,float hab)
{
-x = xx; y = yy;
-pPatch = patch;
-envVal = 1.0; // default - no effect of any gradient
-envDev = eps = 0.0;
-habitats.push_back(hab);
-visits = 0;
-smsData = 0;
+ x = xx; y = yy;
+ pPatch = patch;
+ envVal = 1.0; // default - no effect of any gradient
+ envDev = eps = 0.0;
+ habitats.push_back(hab);
+ visits = 0;
+ smsData = 0;
}
Cell::~Cell() {
-#if RSDEBUG
-//DEBUGLOG << "Cell::~Cell(): this = " << this << " smsData = " << smsData << endl;
-#endif
-habIxx.clear();
-habitats.clear();
-if (smsData != 0) {
- if (smsData->effcosts != 0) delete smsData->effcosts;
- delete smsData;
-}
+ habIxx.clear();
+ habitats.clear();
+ if (smsData != 0) {
+ if (smsData->effcosts != 0) delete smsData->effcosts;
+ delete smsData;
+ }
+demoScalings.clear();
+
#if RSDEBUG
//DEBUGLOG << "Cell::~Cell(): deleted" << endl;
#endif
}
void Cell::setHabIndex(short hx) {
-#if RSDEBUG
-//DebugGUI(("Cell::setHabIndex(): this=" + Int2Str((int)this)
-// + " x=" + Int2Str(x) + " y=" + Int2Str(y)
-// + " habIx=" + Int2Str(habIx)
-//).c_str());
-#endif
-if (hx < 0) habIxx.push_back(0);
-else habIxx.push_back(hx);
+ if (hx < 0) habIxx.push_back(0);
+ else habIxx.push_back(hx);
}
void Cell::changeHabIndex(short ix,short hx) {
-if (ix >= 0 && ix < (short)habIxx.size() && hx >= 0) habIxx[ix] = hx;
-else habIxx[ix] = 0;
+ if (ix >= 0 && ix < (short)habIxx.size() && hx >= 0) habIxx[ix] = hx;
+ else habIxx[ix] = 0;
}
int Cell::getHabIndex(int ix) {
-if (ix < 0 || ix >= (int)habIxx.size())
- // nodata cell OR should not occur, but treat as such
- return -1;
-else return habIxx[ix];
+ if (ix < 0 || ix >= (int)habIxx.size())
+ // nodata cell OR should not occur, but treat as such
+ return -1;
+ else return habIxx[ix];
}
int Cell::nHabitats(void) {
-int nh = (int)habIxx.size();
-if ((int)habitats.size() > nh) nh = (int)habitats.size();
-return nh;
+ int nh = (int)habIxx.size();
+ if ((int)habitats.size() > nh) nh = (int)habitats.size();
+ return nh;
}
void Cell::setHabitat(float q) {
-if (q >= 0.0 && q <= 100.0) habitats.push_back(q);
-else habitats.push_back(0.0);
+ if (q >= 0.0 && q <= 100.0) habitats.push_back(q);
+ else habitats.push_back(0.0);
}
float Cell::getHabitat(int ix) {
-if (ix < 0 || ix >= (int)habitats.size())
- // nodata cell OR should not occur, but treat as such
- return -1.0;
-else return habitats[ix];
+ if (ix < 0 || ix >= (int)habitats.size())
+ // nodata cell OR should not occur, but treat as such
+ return -1.0;
+ else return habitats[ix];
}
-void Cell::setPatch(intptr p) {
-pPatch = p;
+void Cell::setPatch(Patch *p) {
+ pPatch = p;
}
-intptr Cell::getPatch(void)
+Patch *Cell::getPatch(void)
{
-#if RSDEBUG
-//DebugGUI(("Cell::getPatch(): this=" + Int2Str((int)this)
-// + " x=" + Int2Str(x) + " y=" + Int2Str(y)
-// + " habIxx[0]=" + Int2Str(habIxx[0]) + " pPatch=" + Int2Str(pPatch)
-//).c_str());
-#endif
-return pPatch;
+ return pPatch;
}
locn Cell::getLocn(void) { locn q; q.x = x; q.y = y; return q; }
@@ -134,79 +115,104 @@ void Cell::setEnvDev(float d) { envDev = d; }
float Cell::getEnvDev(void) { return envDev; }
void Cell::setEnvVal(float e) {
-if (e >= 0.0) envVal = e;
+ if (e >= 0.0) envVal = e;
}
float Cell::getEnvVal(void) { return envVal; }
void Cell::updateEps(float ac,float randpart) {
-eps = eps*ac + randpart;
+ eps = eps * ac + randpart;
}
float Cell::getEps(void) { return eps; }
// Functions to handle costs for SMS
+#ifdef _OPENMP
+std::unique_lock Cell::lockCost() {
+ return std::unique_lock(cost_mutex);
+}
+#endif
+
int Cell::getCost(void) {
-int c;
-if (smsData == 0) c = 0; // costs not yet set up
-else c = smsData->cost;
-return c;
+ int c;
+ if (smsData == 0) c = 0; // costs not yet set up
+ else c = smsData->cost;
+ return c;
}
void Cell::setCost(int c) {
-if (smsData == 0) {
- smsData = new smscosts;
- smsData->effcosts = 0;
-}
-smsData->cost = c;
+ if (smsData == 0) {
+ smsData = new smscosts;
+ smsData->effcosts = 0;
+ }
+ smsData->cost = c;
}
// Reset the cost and the effective cost of the cell
void Cell::resetCost(void) {
-if (smsData != 0) { resetEffCosts(); delete smsData; }
-smsData = 0;
+ if (smsData != 0) { resetEffCosts(); delete smsData; }
+ smsData = 0;
}
array3x3f Cell::getEffCosts(void) {
-array3x3f a;
-if (smsData == 0 || smsData->effcosts == 0) { // effective costs have not been calculated
- for (int i = 0; i < 3; i++) {
- for (int j = 0; j < 3; j++) {
- a.cell[i][j] = -1.0;
+ array3x3f a;
+ if (smsData == 0 || smsData->effcosts == 0) { // effective costs have not been calculated
+ for (int i = 0; i < 3; i++) {
+ for (int j = 0; j < 3; j++) {
+ a.cell[i][j] = -1.0;
+ }
}
}
-}
-else
- a = *smsData->effcosts;
-return a;
+ else
+ a = *smsData->effcosts;
+ return a;
}
-void Cell::setEffCosts(array3x3f a) {
-if (smsData->effcosts == 0) smsData->effcosts = new array3x3f;
-*smsData->effcosts = a;
+void Cell::setEffCosts(array3x3f a) {
+ if (smsData->effcosts == 0) smsData->effcosts = new array3x3f;
+ *smsData->effcosts = a;
}
// Reset the effective cost, but not the cost, of the cell
void Cell::resetEffCosts(void) {
-if (smsData != 0) {
- if (smsData->effcosts != 0) {
- delete smsData->effcosts;
- smsData->effcosts = 0;
+ if (smsData != 0) {
+ if (smsData->effcosts != 0) {
+ delete smsData->effcosts;
+ smsData->effcosts = 0;
+ }
}
}
-}
void Cell::resetVisits(void) { visits = 0; }
void Cell::incrVisits(void) { visits++; }
unsigned long int Cell::getVisits(void) { return visits; }
+
+void Cell::addchgDemoScaling(std::vector ds) {
+ std::for_each(ds.begin(), ds.end(), [](float& perc){ if(perc < 0.0 || perc > 100.0) perc=100; });
+ demoScalings.push_back(ds);
+ return;
+}
+
+std::vector Cell::getDemoScaling(short chgyear) {
+ if (chgyear < 0 || chgyear >= (int)demoScalings.size()) {
+ std::vector ret(1, -1);
+ return ret;
+ }
+ else return demoScalings[chgyear];
+}
+
+
+
//---------------------------------------------------------------------------
// Initial species distribution cell functions
DistCell::DistCell(int xx,int yy) {
-x = xx; y = yy; initialise = false;
+ x = xx;
+ y = yy;
+ initialise = false;
}
DistCell::~DistCell() {
@@ -214,18 +220,21 @@ DistCell::~DistCell() {
}
void DistCell::setCell(bool init) {
-initialise = init;
+ initialise = init;
}
bool DistCell::toInitialise(locn loc) {
-if (loc.x == x && loc.y == y) return initialise;
-else return false;
+ if (loc.x == x && loc.y == y) return initialise;
+ else return false;
}
bool DistCell::selected(void) { return initialise; }
locn DistCell::getLocn(void) {
-locn loc; loc.x = x; loc.y = y; return loc;
+ locn loc;
+ loc.x = x;
+ loc.y = y;
+ return loc;
}
//---------------------------------------------------------------------------
diff --git a/Cell.h b/Cell.h
index 5382a1e..2dc28f8 100644
--- a/Cell.h
+++ b/Cell.h
@@ -1,74 +1,84 @@
/*----------------------------------------------------------------------------
- *
- * Copyright (C) 2020 Greta Bocedi, Stephen C.F. Palmer, Justin M.J. Travis, Anne-Kathleen Malchow, Damaris Zurell
- *
+ *
+ * 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 .
- *
+ *
--------------------------------------------------------------------------*/
-
-
-/*------------------------------------------------------------------------------
-RangeShifter v2.0 Cell
-Implements the following classes:
+ /*------------------------------------------------------------------------------
-Cell - Landscape cell
+ RangeShifter v2.0 Cell
-DistCell - Initial species distribution cell
+ Implements the following classes:
-For full details of RangeShifter, please see:
-Bocedi G., Palmer S.C.F., Peer G., Heikkinen R.K., Matsinos Y.G., Watts K.
-and Travis J.M.J. (2014). RangeShifter: a platform for modelling spatial
-eco-evolutionary dynamics and species responses to environmental changes.
-Methods in Ecology and Evolution, 5, 388-396. doi: 10.1111/2041-210X.12162
+ Cell - Landscape cell
-Authors: Greta Bocedi & Steve Palmer, University of Aberdeen
+ DistCell - Initial species distribution cell
-Last updated: 14 January 2021 by Steve Palmer
+ For full details of RangeShifter, please see:
+ Bocedi G., Palmer S.C.F., Pe’er G., Heikkinen R.K., Matsinos Y.G., Watts K.
+ and Travis J.M.J. (2014). RangeShifter: a platform for modelling spatial
+ eco-evolutionary dynamics and species’ responses to environmental changes.
+ Methods in Ecology and Evolution, 5, 388-396. doi: 10.1111/2041-210X.12162
-------------------------------------------------------------------------------*/
+ Authors: Greta Bocedi & Steve Palmer, University of Aberdeen
+
+ Last updated: 14 January 2021 by Steve Palmer
+
+ ------------------------------------------------------------------------------*/
#ifndef CellH
#define CellH
+
+#include
+
#include
using namespace std;
#include "Parameters.h"
+#ifdef _OPENMP
+#include
+#include
+#endif
+
//---------------------------------------------------------------------------
+class Patch; // Forward-declaration of the Patch class
+
struct array3x3f { float cell[3][3]; }; // neighbourhood cell array (SMS)
-struct smscosts { int cost; array3x3f *effcosts; }; // cell costs for SMS
+struct smscosts { int cost; array3x3f* effcosts; }; // cell costs for SMS
// Landscape cell
-class Cell{
+class Cell {
public:
Cell( // Constructor for habitat codes
int, // x co-ordinate
int, // y co-ordinate
- intptr, // pointer (cast as integer) to the Patch to which Cell belongs
+ Patch *, // pointer to the Patch to which Cell belongs
int // habitat index number
);
Cell( // Constructor for habitat % cover or habitat quality
int, // x co-ordinate
int, // y co-ordinate
- intptr, // pointer (cast as integer) to the Patch to which Cell belongs
+ Patch *, // pointer to the Patch to which Cell belongs
float // habitat proportion or cell quality score
);
~Cell();
@@ -90,9 +100,9 @@ class Cell{
int // habitat index number / landscape change number
);
void setPatch(
- intptr // pointer (cast as integer) to the Patch to which Cell belongs
+ Patch * // pointer to the Patch to which Cell belongs
);
- intptr getPatch(void);
+ Patch *getPatch(void);
locn getLocn(void);
void setEnvDev(
float // local environmental deviation
@@ -110,6 +120,9 @@ class Cell{
void setCost(
int // cost value for SMS
);
+#ifdef _OPENMP
+ std::unique_lock lockCost(void);
+#endif
int getCost(void);
void resetCost(void);
array3x3f getEffCosts(void);
@@ -121,29 +134,44 @@ class Cell{
void incrVisits(void);
unsigned long int getVisits(void);
+ void addchgDemoScaling(std::vector);
+ void setDemoScaling(std::vector, short);
+ std::vector getDemoScaling(short);
+
+
private:
- int x,y; // cell co-ordinates
- intptr pPatch; // pointer (cast as integer) to the Patch to which cell belongs
+ int x, y; // cell co-ordinates
+ Patch *pPatch; // pointer to the Patch to which cell belongs
// NOTE: THE FOLLOWING ENVIRONMENTAL VARIABLES COULD BE COMBINED IN A STRUCTURE
// AND ACCESSED BY A POINTER ...
float envVal; // environmental value, representing one of:
- // gradient in K, r or extinction probability
+ // gradient in K, r or extinction probability
float envDev; // local environmental deviation (static, in range -1.0 to +1.0)
float eps; // local environmental stochasticity (epsilon) (dynamic, from N(0,std))
+#ifdef _OPENMP
+ std::atomic visits; // no. of times square is visited by dispersers
+#else
unsigned long int visits; // no. of times square is visited by dispersers
- smscosts *smsData;
+#endif
+ smscosts* smsData;
vector habIxx; // habitat indices (rasterType=0)
- // NB initially, habitat codes are loaded, then converted to index nos.
- // once landscape is fully loaded
+ // NB initially, habitat codes are loaded, then converted to index nos.
+ // once landscape is fully loaded
vector habitats; // habitat proportions (rasterType=1) or quality (rasterType=2)
+
+ std::vector> demoScalings; // demographic scaling layers (only if rasterType==2)
+
+#ifdef _OPENMP
+ std::mutex cost_mutex;
+#endif
};
//---------------------------------------------------------------------------
// Initial species distribution cell
-class DistCell{
+class DistCell {
public:
DistCell(
int, // x co-ordinate
@@ -160,15 +188,11 @@ class DistCell{
locn getLocn(void);
private:
- int x,y; // cell co-ordinates
+ int x, y; // cell co-ordinates
bool initialise; // cell is to be initialised
};
-#if RSDEBUG
-extern void DebugGUI(string);
-#endif
-
//---------------------------------------------------------------------------
#endif
diff --git a/Community.cpp b/Community.cpp
index 286dc6a..46c8970 100644
--- a/Community.cpp
+++ b/Community.cpp
@@ -1,6 +1,7 @@
+#include "Community.h"
/*----------------------------------------------------------------------------
*
- * Copyright (C) 2020 Greta Bocedi, Stephen C.F. Palmer, Justin M.J. Travis, Anne-Kathleen Malchow, Damaris Zurell
+ * 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.
*
@@ -24,18 +25,41 @@
#include "Community.h"
+#ifdef _OPENMP
+#ifdef __has_include
+#if __has_include()
+#include
+#endif
+#endif
+#include
+#if __cpp_lib_barrier >= 201907L && __cpp_lib_optional >= 201606L
+#define HAS_BARRIER_LIB
+#include
+#include
+#else
+#include
+#include
+#endif
+#include
+#endif // _OPENMP
+
//---------------------------------------------------------------------------
ofstream outrange;
ofstream outoccup, outsuit;
ofstream outtraitsrows;
+ofstream ofsGenes;
+ofstream outwcfstat;
+ofstream outperlocusfstat;
+ofstream outpairwisefst;
//---------------------------------------------------------------------------
Community::Community(Landscape* pLand) {
pLandscape = pLand;
indIx = 0;
+ pNeutralStatistics = 0;
}
Community::~Community(void) {
@@ -54,13 +78,11 @@ SubCommunity* Community::addSubComm(Patch* pPch, int num) {
void Community::initialise(Species* pSpecies, int year)
{
-
int nsubcomms, npatches, ndistcells, spratio, patchnum, rr = 0;
locn distloc;
patchData pch;
patchLimits limits;
- intptr ppatch, subcomm;
- std::vector subcomms;
+ std::vector subcomms;
std::vector selected;
SubCommunity* pSubComm;
Patch* pPatch;
@@ -72,16 +94,6 @@ void Community::initialise(Species* pSpecies, int year)
spratio = ppLand.spResol / ppLand.resol;
-#if RSDEBUG
- DEBUGLOG << endl << "Community::initialise(): this=" << this
- << " seedType=" << init.seedType << " freeType=" << init.freeType
- << " minSeedX=" << init.minSeedX << " minSeedY=" << init.minSeedY
- << " maxSeedX=" << init.maxSeedX << " maxSeedY=" << init.maxSeedY
- << " indsFile=" << init.indsFile
- << " nsubcomms=" << nsubcomms << " spratio=" << spratio
- << endl;
-#endif
-
switch (init.seedType) {
case 0: // free initialisation
@@ -138,7 +150,7 @@ void Community::initialise(Species* pSpecies, int year)
}
for (int i = 0; i < npatches; i++) {
if (selected[i]) {
- pSubComm = (SubCommunity*)subcomms[i];
+ pSubComm = subcomms[i];
pSubComm->setInitial(true);
}
}
@@ -155,14 +167,11 @@ void Community::initialise(Species* pSpecies, int year)
if (patchnum != 0) {
if (pch.pPatch->getK() > 0.0)
{ // patch is suitable
- subcomm = pch.pPatch->getSubComm();
- if (subcomm == 0) {
+ pSubComm = pch.pPatch->getSubComm();
+ if (pSubComm == nullptr) {
// create a sub-community in the patch
pSubComm = addSubComm(pch.pPatch, patchnum);
}
- else {
- pSubComm = (SubCommunity*)subcomm;
- }
pSubComm->setInitial(true);
}
}
@@ -211,13 +220,11 @@ void Community::initialise(Species* pSpecies, int year)
for (int y = 0; y < spratio; y++) {
pCell = pLandscape->findCell(distloc.x * spratio + x, distloc.y * spratio + y);
if (pCell != 0) { // not a no-data cell
- ppatch = pCell->getPatch();
- if (ppatch != 0) {
- pPatch = (Patch*)ppatch;
+ pPatch = pCell->getPatch();
+ if (pPatch != nullptr) {
if (pPatch->getSeqNum() != 0) { // not the matrix patch
- subcomm = pPatch->getSubComm();
- if (subcomm != 0) {
- pSubComm = (SubCommunity*)subcomm;
+ pSubComm = pPatch->getSubComm();
+ if (pSubComm != nullptr) {
pSubComm->setInitial(true);
}
}
@@ -246,7 +253,8 @@ void Community::initialise(Species* pSpecies, int year)
indIx = 0; // reset index for initial individuals
}
else { // add any initial individuals for the current year
- initInd iind; iind.year = 0;
+ initInd iind = initInd();
+ iind.year = 0;
int ninds = paramsInit->numInitInds();
while (indIx < ninds && iind.year <= year) {
iind = paramsInit->getInitInd(indIx);
@@ -256,14 +264,11 @@ void Community::initialise(Species* pSpecies, int year)
pPatch = pLandscape->findPatch(iind.patchID);
if (pPatch->getK() > 0.0)
{ // patch is suitable
- subcomm = pPatch->getSubComm();
- if (subcomm == 0) {
+ pSubComm = pPatch->getSubComm();
+ if (pSubComm == nullptr) {
// create a sub-community in the patch
pSubComm = addSubComm(pPatch, iind.patchID);
}
- else {
- pSubComm = (SubCommunity*)subcomm;
- }
pSubComm->initialInd(pLandscape, pSpecies, pPatch, pPatch->getRandomCell(), indIx);
}
}
@@ -271,19 +276,15 @@ void Community::initialise(Species* pSpecies, int year)
else { // cell-based model
pCell = pLandscape->findCell(iind.x, iind.y);
if (pCell != 0) {
- intptr ppatch = pCell->getPatch();
- if (ppatch != 0) {
- pPatch = (Patch*)ppatch;
+ pPatch = pCell->getPatch();
+ if (pPatch != nullptr) {
if (pPatch->getK() > 0.0)
{ // patch is suitable
- subcomm = pPatch->getSubComm();
- if (subcomm == 0) {
+ pSubComm = pPatch->getSubComm();
+ if (pSubComm == nullptr) {
// create a sub-community in the patch
pSubComm = addSubComm(pPatch, iind.patchID);
}
- else {
- pSubComm = (SubCommunity*)subcomm;
- }
pSubComm->initialInd(pLandscape, pSpecies, pPatch, pCell, indIx);
}
}
@@ -308,18 +309,11 @@ void Community::initialise(Species* pSpecies, int year)
} // end of switch (init.seedType)
-#if RSDEBUG
- DEBUGLOG << "Community::initialise(): this=" << this
- << " nsubcomms=" << nsubcomms
- << endl;
-#endif
-
}
// Add manually selected patches/cells to the selected set for initialisation
void Community::addManuallySelected(void) {
int npatches;
- intptr subcomm, patch;
locn initloc;
Cell* pCell;
Patch* pPatch;
@@ -328,19 +322,14 @@ void Community::addManuallySelected(void) {
landParams ppLand = pLandscape->getLandParams();
npatches = pLandscape->initCellCount(); // no. of patches/cells specified
-#if RSDEBUG
- DEBUGLOG << "Community::addManuallySelected(): this = " << this
- << " npatches = " << npatches << endl;
-#endif
// identify sub-communities to be initialised
if (ppLand.patchModel) {
for (int i = 0; i < npatches; i++) {
initloc = pLandscape->getInitCell(i); // patch number held in x-coord of list
pPatch = pLandscape->findPatch(initloc.x);
if (pPatch != 0) {
- subcomm = pPatch->getSubComm();
- if (subcomm != 0) {
- pSubComm = (SubCommunity*)subcomm;
+ pSubComm = pPatch->getSubComm();
+ if (pSubComm != nullptr) {
pSubComm->setInitial(true);
}
}
@@ -353,29 +342,11 @@ void Community::addManuallySelected(void) {
&& initloc.y >= 0 && initloc.y < ppLand.dimY) {
pCell = pLandscape->findCell(initloc.x, initloc.y);
if (pCell != 0) { // not no-data cell
- patch = pCell->getPatch();
-#if RSDEBUG
- DEBUGLOG << "Community::initialise(): i = " << i
- << " x = " << initloc.x << " y = " << initloc.y
- << " pCell = " << pCell << " patch = " << patch
- << endl;
-#endif
- if (patch != 0) {
- pPatch = (Patch*)patch;
- subcomm = pPatch->getSubComm();
-#if RSDEBUG
- DEBUGLOG << "Community::initialise(): i = " << i
- << " pPatch = " << pPatch << " subcomm = " << subcomm
- << endl;
-#endif
- if (subcomm != 0) {
- pSubComm = (SubCommunity*)subcomm;
+ pPatch = pCell->getPatch();
+ if (pPatch != nullptr) {
+ pSubComm = pPatch->getSubComm();
+ if (pSubComm != nullptr) {
pSubComm->setInitial(true);
-#if RSDEBUG
- DEBUGLOG << "Community::initialise(): i = " << i
- << " pSubComm = " << pSubComm
- << endl;
-#endif
}
}
}
@@ -386,6 +357,7 @@ void Community::addManuallySelected(void) {
void Community::resetPopns(void) {
int nsubcomms = (int)subComms.size();
+ #pragma omp parallel for schedule(static,128)
for (int i = 0; i < nsubcomms; i++) { // all sub-communities
subComms[i]->resetPopns();
}
@@ -413,49 +385,93 @@ void Community::patchChanges(void) {
void Community::reproduction(int yr)
{
- float eps = 0.0; // epsilon for environmental stochasticity
landParams land = pLandscape->getLandParams();
envStochParams env = paramsStoch->getStoch();
+ float eps = 0.0; // epsilon for environmental stochasticity
+ if (env.stoch && !env.local) { // global stochasticty
+ eps = pLandscape->getGlobalStoch(yr);
+ }
int nsubcomms = (int)subComms.size();
-#if RSDEBUG
- DEBUGLOG << "Community::reproduction(): this=" << this
- << " nsubcomms=" << nsubcomms << endl;
-#endif
+ #pragma omp parallel for schedule(static,128)
for (int i = 0; i < nsubcomms; i++) { // all sub-communities
- if (env.stoch) {
- if (!env.local) { // global stochasticty
- eps = pLandscape->getGlobalStoch(yr);
- }
- }
subComms[i]->reproduction(land.resol, eps, land.rasterType, land.patchModel);
}
-#if RSDEBUG
- DEBUGLOG << "Community::reproduction(): finished" << endl;
-#endif
}
void Community::emigration(void)
{
- int nsubcomms = (int)subComms.size();
-#if RSDEBUG
- DEBUGLOG << "Community::emigration(): this=" << this
- << " nsubcomms=" << nsubcomms << endl;
-#endif
+ int nsubcomms = static_cast(subComms.size());
+ #pragma omp parallel for schedule(static, 128)
for (int i = 0; i < nsubcomms; i++) { // all sub-communities
subComms[i]->emigration();
}
-#if RSDEBUG
- DEBUGLOG << "Community::emigration(): finished" << endl;
-#endif
}
-#if RS_RCPP // included also SEASONAL
+#ifdef _OPENMP
+#ifdef HAS_BARRIER_LIB
+typedef std::optional> split_barrier;
+#else
+class split_barrier {
+private:
+ std::mutex m;
+ std::condition_variable cv;
+ int threads_in_section;
+ int total_threads;
+ bool may_enter;
+ bool may_leave;
+
+public:
+ split_barrier():
+ threads_in_section(0),
+ may_enter(false),
+ may_leave(false)
+ {}
+
+ void emplace(int threads) {
+ std::lock_guard lock(m);
+ total_threads = threads;
+ may_enter = true;
+ }
+
+ void enter() {
+ std::unique_lock lock(m);
+ cv.wait(lock, [this]{return may_enter;});
+ if (++threads_in_section == total_threads) {
+ may_enter = false;
+ may_leave = true;
+ lock.unlock();
+ cv.notify_all();
+ }
+ }
+
+ void leave() {
+ std::unique_lock lock(m);
+ cv.wait(lock, [this]{return may_leave;});
+ if (--threads_in_section == 0) {
+ may_leave = false;
+ may_enter = true;
+ lock.unlock();
+ cv.notify_all();
+ }
+ }
+};
+#endif // HAS_BARRIER_LIB
+#endif // _OPENMP
+
+#if RS_RCPP
void Community::dispersal(short landIx, short nextseason)
#else
void Community::dispersal(short landIx)
-#endif // SEASONAL || RS_RCPP
+#endif // RS_RCPP
{
+#ifdef _OPENMP
+ std::atomic nbStillDispersing;
+ split_barrier barrier;
+#else
+ int nbStillDispersing;
+#endif // _OPENMP
+
#if RSDEBUG
int t0, t1, t2;
t0 = time(0);
@@ -466,43 +482,89 @@ void Community::dispersal(short landIx)
int nsubcomms = (int)subComms.size();
// initiate dispersal - all emigrants leave their natal community and join matrix community
SubCommunity* matrix = subComms[0]; // matrix community is always the first
- for (int i = 0; i < nsubcomms; i++) { // all populations
- subComms[i]->initiateDispersal(matrix);
- }
-#if RSDEBUG
- t1 = time(0);
- DEBUGLOG << "Community::dispersal(): this=" << this
- << " nsubcomms=" << nsubcomms << " initiation time=" << t1 - t0 << endl;
-#endif
+#pragma omp parallel
+ {
+ std::map> disperserPool;
- // dispersal is undertaken by all individuals now in the matrix patch
- // (even if not physically in the matrix)
- int ndispersers = 0;
- do {
- for (int i = 0; i < nsubcomms; i++) { // all populations
- subComms[i]->resetPossSettlers();
+ // All individuals in the matrix disperse again
+ // (= unsettled dispersers from previous generation)
+ matrix->disperseMatrix(disperserPool);
+
+ // Recruit new emigrants
+#pragma omp for schedule(static,128) nowait
+ for (int i = 0; i < nsubcomms; i++) {
+ subComms[i]->recruitDispersers(disperserPool);
}
-#if RS_RCPP // included also SEASONAL
- ndispersers = matrix->transfer(pLandscape, landIx, nextseason);
+
+#ifdef _OPENMP
+#pragma omp single
+ barrier.emplace(omp_get_num_threads());
+#endif // _OPENMP
+
+
+ //
+ do {
+#pragma omp for schedule(guided)
+ for (int i = 0; i < nsubcomms; i++) {
+ subComms[i]->resetPossSettlers();
+ }
+ int localNbDispersers = matrix->resolveTransfer(disperserPool, pLandscape, landIx);
+#pragma omp single nowait
+ nbStillDispersing = 0;
+#pragma omp barrier
+#if RS_RCPP
+ localNbDispersers += matrix->resolveSettlement(disperserPool, pLandscape, nextseason);
#else
- ndispersers = matrix->transfer(pLandscape, landIx);
-#endif // SEASONAL || RS_RCPP
- matrix->completeDispersal(pLandscape, sim.outConnect);
- } while (ndispersers > 0);
+ localNbDispersers += matrix->resolveSettlement(disperserPool, pLandscape);
+#endif // RS_RCPP
+ nbStillDispersing += localNbDispersers;
+
+#ifdef _OPENMP
+#ifdef HAS_BARRIER_LIB
+ std::barrier<>::arrival_token token = barrier->arrive();
+#else
+ barrier.enter();
+#endif // HAS_BARRIER_LIB
+#endif // _OPENMP
+
+ matrix->completeDispersal(disperserPool, pLandscape, sim.outConnect);
+
+#ifdef _OPENMP
+#ifdef HAS_BARRIER_LIB
+ barrier->wait(std::move(token));
+#else
+ barrier.leave();
+#endif // HAS_BARRIER_LIB
+#endif // _OPENMP
+
+ } while (nbStillDispersing > 0);
+
+ // All unsettled dispersers are stored in matrix until next generation
+ for (auto & it : disperserPool) {
+ Species* const& pSpecies = it.first;
+ vector& inds = it.second;
+ matrix->recruitMany(inds, pSpecies);
+ }
+ }
-#if RSDEBUG
- DEBUGLOG << "Community::dispersal(): matrix=" << matrix << endl;
- t2 = time(0);
- DEBUGLOG << "Community::dispersal(): transfer time=" << t2 - t1 << endl;
-#endif
}
-void Community::survival(short part, short option0, short option1)
+void Community::survival0(short option0, short option1)
{
int nsubcomms = (int)subComms.size();
+ #pragma omp parallel for schedule(static,128)
for (int i = 0; i < nsubcomms; i++) { // all communities (including in matrix)
- subComms[i]->survival(part, option0, option1);
+ subComms[i]->survival0(option0, option1);
+ }
+}
+
+void Community::survival1()
+{
+ int nsubcomms = (int)subComms.size();
+ #pragma omp parallel for schedule(static,128)
+ for (int i = 0; i < nsubcomms; i++) { // all communities (including in matrix)
+ subComms[i]->survival1();
}
}
@@ -553,17 +615,12 @@ void Community::createOccupancy(int nrows, int reps) {
void Community::updateOccupancy(int row, int rep)
{
-#if RSDEBUG
- DEBUGLOG << "Community::updateOccupancy(): row=" << row << endl;
-#endif
int nsubcomms = (int)subComms.size();
for (int i = 0; i < nsubcomms; i++) {
subComms[i]->updateOccupancy(row);
}
-
commStats s = getStats();
occSuit[row][rep] = (float)s.occupied / (float)s.suitable;
-
}
void Community::deleteOccupancy(int nrows) {
@@ -583,7 +640,7 @@ void Community::deleteOccupancy(int nrows) {
// Determine range margins
commStats Community::getStats(void)
{
- commStats s;
+ commStats s = commStats();
landParams ppLand = pLandscape->getLandParams();
s.ninds = s.nnonjuvs = s.suitable = s.occupied = 0;
s.minX = ppLand.maxX; s.minY = ppLand.maxY; s.maxX = s.maxY = 0;
@@ -616,9 +673,14 @@ commStats Community::getStats(void)
// Functions to control production of output files
+// Close population file
+bool Community::outPopFinishLandscape() {
+ return subComms[0]->outPopFinishLandscape();
+}
+
// Open population file and write header record
-bool Community::outPopHeaders(Species* pSpecies, int option) {
- return subComms[0]->outPopHeaders(pLandscape, pSpecies, option);
+bool Community::outPopStartLandscape(Species* pSpecies) {
+ return subComms[0]->outPopStartLandscape(pLandscape, pSpecies);
}
// Write records to population file
@@ -632,53 +694,36 @@ void Community::outPop(int rep, int yr, int gen)
}
+// Close individuals file
+void Community::outIndsFinishReplicate() {
+ subComms[0]->outIndsFinishReplicate();
+}
-// Write records to individuals file
-void Community::outInds(int rep, int yr, int gen, int landNr) {
+// Open individuals file and write header record
+void Community::outIndsStartReplicate(int rep, int landNr) {
+ subComms[0]->outIndsStartReplicate(pLandscape, rep, landNr);
+}
- if (landNr >= 0) { // open the file
- subComms[0]->outInds(pLandscape, rep, yr, gen, landNr);
- return;
- }
- if (landNr == -999) { // close the file
- subComms[0]->outInds(pLandscape, rep, yr, gen, -999);
- return;
- }
+// Write records to individuals file
+void Community::outIndividuals(int rep, int yr, int gen) {
// generate output for each sub-community (patch) in the community
int nsubcomms = (int)subComms.size();
for (int i = 0; i < nsubcomms; i++) { // all sub-communities
- subComms[i]->outInds(pLandscape, rep, yr, gen, landNr);
+ subComms[i]->outIndividuals(pLandscape, rep, yr, gen);
}
}
-// Write records to genetics file
-void Community::outGenetics(int rep, int yr, int gen, int landNr) {
- //landParams ppLand = pLandscape->getLandParams();
- if (landNr >= 0) { // open the file
- subComms[0]->outGenetics(rep, yr, gen, landNr);
- return;
- }
- if (landNr == -999) { // close the file
- subComms[0]->outGenetics(rep, yr, gen, landNr);
- return;
- }
- // generate output for each sub-community (patch) in the community
- int nsubcomms = (int)subComms.size();
- for (int i = 0; i < nsubcomms; i++) { // all sub-communities
- subComms[i]->outGenetics(rep, yr, gen, landNr);
- }
+// Close range file
+bool Community::outRangeFinishLandscape()
+{
+ if (outrange.is_open()) outrange.close();
+ outrange.clear();
+ return true;
}
// Open range file and write header record
-bool Community::outRangeHeaders(Species* pSpecies, int landNr)
+bool Community::outRangeStartLandscape(Species* pSpecies, int landNr)
{
-
- if (landNr == -999) { // close the file
- if (outrange.is_open()) outrange.close();
- outrange.clear();
- return true;
- }
-
string name;
landParams ppLand = pLandscape->getLandParams();
envStochParams env = paramsStoch->getStoch();
@@ -686,33 +731,21 @@ bool Community::outRangeHeaders(Species* pSpecies, int landNr)
// NEED TO REPLACE CONDITIONAL COLUMNS BASED ON ATTRIBUTES OF ONE SPECIES TO COVER
// ATTRIBUTES OF *ALL* SPECIES AS DETECTED AT MODEL LEVEL
- demogrParams dem = pSpecies->getDemogr();
- stageParams sstruct = pSpecies->getStage();
- emigRules emig = pSpecies->getEmig();
- trfrRules trfr = pSpecies->getTrfr();
+ demogrParams dem = pSpecies->getDemogrParams();
+ stageParams sstruct = pSpecies->getStageParams();
+ emigRules emig = pSpecies->getEmigRules();
+ transferRules trfr = pSpecies->getTransferRules();
settleType sett = pSpecies->getSettle();
-#if RSDEBUG
- DEBUGLOG << "Community::outRangeHeaders(): simulation=" << sim.simulation
- << " sim.batchMode=" << sim.batchMode
- << " landNr=" << landNr << endl;
-#endif
-
if (sim.batchMode) {
name = paramsSim->getDir(2)
-#if RS_RCPP
- + "Batch" + Int2Str(sim.batchNum) + "_"
- + "Sim" + Int2Str(sim.simulation) + "_Land"
- + Int2Str(landNr)
-#else
- + "Batch" + Int2Str(sim.batchNum) + "_"
- + "Sim" + Int2Str(sim.simulation) + "_Land"
- + Int2Str(landNr)
-#endif
+ + "Batch" + to_string(sim.batchNum) + "_"
+ + "Sim" + to_string(sim.simulation) + "_Land"
+ + to_string(landNr)
+ "_Range.txt";
}
else {
- name = paramsSim->getDir(2) + "Sim" + Int2Str(sim.simulation) + "_Range.txt";
+ name = paramsSim->getDir(2) + "Sim" + to_string(sim.simulation) + "_Range.txt";
}
outrange.open(name.c_str());
outrange << "Rep\tYear\tRepSeason";
@@ -749,7 +782,7 @@ bool Community::outRangeHeaders(Species* pSpecies, int landNr)
}
}
if (trfr.indVar) {
- if (trfr.moveModel) {
+ if (trfr.usesMovtProc) {
if (trfr.moveType == 1) {
outrange << "\tmeanDP\tstdDP\tmeanGB\tstdGB";
outrange << "\tmeanAlphaDB\tstdAlphaDB\tmeanBetaDB\tstdBetaDB";
@@ -787,31 +820,21 @@ bool Community::outRangeHeaders(Species* pSpecies, int landNr)
}
}
outrange << endl;
-
-#if RSDEBUG
- DEBUGLOG << "Community::outRangeHeaders(): finished" << endl;
-#endif
-
return outrange.is_open();
}
// Write record to range file
void Community::outRange(Species* pSpecies, int rep, int yr, int gen)
{
-#if RSDEBUG
- DEBUGLOG << "Community::outRange(): rep=" << rep
- << " yr=" << yr << " gen=" << gen << endl;
-#endif
-
landParams ppLand = pLandscape->getLandParams();
envStochParams env = paramsStoch->getStoch();
// NEED TO REPLACE CONDITIONAL COLUMNS BASED ON ATTRIBUTES OF ONE SPECIES TO COVER
// ATTRIBUTES OF *ALL* SPECIES AS DETECTED AT MODEL LEVEL
- demogrParams dem = pSpecies->getDemogr();
- stageParams sstruct = pSpecies->getStage();
- emigRules emig = pSpecies->getEmig();
- trfrRules trfr = pSpecies->getTrfr();
+ demogrParams dem = pSpecies->getDemogrParams();
+ stageParams sstruct = pSpecies->getStageParams();
+ emigRules emig = pSpecies->getEmigRules();
+ transferRules trfr = pSpecies->getTransferRules();
settleType sett = pSpecies->getSettle();
outrange << rep << "\t" << yr << "\t" << gen;
@@ -828,14 +851,14 @@ void Community::outRange(Species* pSpecies, int rep, int yr, int gen)
for (int stg = 1; stg < sstruct.nStages; stg++) {
stagepop = 0;
for (int i = 0; i < nsubcomms; i++) { // all sub-communities
- stagepop += subComms[i]->stagePop(stg);
+ stagepop += subComms[i]->getNbInds(stg);
}
outrange << "\t" << stagepop;
}
// juveniles born in current reproductive season
stagepop = 0;
for (int i = 0; i < nsubcomms; i++) { // all sub-communities
- stagepop += subComms[i]->stagePop(0);
+ stagepop += subComms[i]->getNbInds(0);
}
outrange << "\t" << stagepop;
}
@@ -858,14 +881,11 @@ void Community::outRange(Species* pSpecies, int rep, int yr, int gen)
outrange << "\t0\t0\t0\t0";
if (emig.indVar || trfr.indVar || sett.indVar) { // output trait means
- traitsums ts;
+ traitsums ts = traitsums();
traitsums scts; // sub-community traits
- traitCanvas tcanv;
int ngenes, popsize;
- tcanv.pcanvas[0] = NULL;
-
- for (int i = 0; i < NSEXES; i++) {
+ 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;
@@ -882,8 +902,8 @@ void Community::outRange(Species* pSpecies, int rep, int yr, int gen)
int nsubcomms = (int)subComms.size();
for (int i = 0; i < nsubcomms; i++) { // all sub-communities (incl. matrix)
- scts = subComms[i]->outTraits(tcanv, pLandscape, rep, yr, gen, true);
- for (int j = 0; j < NSEXES; j++) {
+ scts = subComms[i]->outTraits(pLandscape, rep, yr, gen, true);
+ for (int j = 0; j < gMaxNbSexes; j++) {
ts.ninds[j] += scts.ninds[j];
ts.sumD0[j] += scts.sumD0[j]; ts.ssqD0[j] += scts.ssqD0[j];
ts.sumAlpha[j] += scts.sumAlpha[j]; ts.ssqAlpha[j] += scts.ssqAlpha[j];
@@ -959,7 +979,7 @@ void Community::outRange(Species* pSpecies, int rep, int yr, int gen)
}
if (trfr.indVar) {
- if (trfr.moveModel) {
+ if (trfr.usesMovtProc) {
// CURRENTLY INDIVIDUAL VARIATION CANNOT BE SEX-DEPENDENT
ngenes = 1;
}
@@ -1016,7 +1036,7 @@ void Community::outRange(Species* pSpecies, int rep, int yr, int gen)
}
}
}
- if (trfr.moveModel) {
+ if (trfr.usesMovtProc) {
if (trfr.moveType == 1) {
outrange << "\t" << mnDP[0] << "\t" << sdDP[0];
outrange << "\t" << mnGB[0] << "\t" << sdGB[0];
@@ -1108,16 +1128,18 @@ void Community::outRange(Species* pSpecies, int rep, int yr, int gen)
outrange << endl;
}
-// Open occupancy file, write header record and set up occupancy array
-bool Community::outOccupancyHeaders(int option)
+// Close occupancy file
+bool Community::outOccupancyFinishLandscape()
{
- if (option == -999) { // close the files
- if (outsuit.is_open()) outsuit.close();
- if (outoccup.is_open()) outoccup.close();
- outsuit.clear(); outoccup.clear();
- return true;
- }
+ if (outsuit.is_open()) outsuit.close();
+ if (outoccup.is_open()) outoccup.close();
+ outsuit.clear(); outoccup.clear();
+ return true;
+}
+// Open occupancy file, write header record and set up occupancy array
+bool Community::outOccupancyStartLandscape()
+{
string name, nameI;
simParams sim = paramsSim->getSim();
landParams ppLand = pLandscape->getLandParams();
@@ -1125,22 +1147,22 @@ bool Community::outOccupancyHeaders(int option)
name = paramsSim->getDir(2);
if (sim.batchMode) {
- name += "Batch" + Int2Str(sim.batchNum) + "_";
- name += "Sim" + Int2Str(sim.simulation) + "_Land" + Int2Str(ppLand.landNum);
+ name += "Batch" + to_string(sim.batchNum) + "_";
+ name += "Sim" + to_string(sim.simulation) + "_Land" + to_string(ppLand.landNum);
}
else
- name += "Sim" + Int2Str(sim.simulation);
+ name += "Sim" + to_string(sim.simulation);
name += "_Occupancy_Stats.txt";
outsuit.open(name.c_str());
outsuit << "Year\tMean_OccupSuit\tStd_error" << endl;
name = paramsSim->getDir(2);
if (sim.batchMode) {
- name += "Batch" + Int2Str(sim.batchNum) + "_";
- name += "Sim" + Int2Str(sim.simulation) + "_Land" + Int2Str(ppLand.landNum);
+ name += "Batch" + to_string(sim.batchNum) + "_";
+ name += "Sim" + to_string(sim.simulation) + "_Land" + to_string(ppLand.landNum);
}
else
- name += "Sim" + Int2Str(sim.simulation);
+ name += "Sim" + to_string(sim.simulation);
name += "_Occupancy.txt";
outoccup.open(name.c_str());
if (ppLand.patchModel) {
@@ -1181,9 +1203,10 @@ void Community::outOccupancy(void) {
}
}
-void Community::outOccSuit(bool view) {
+void Community::outOccSuit() {
double sum, ss, mean, sd, se;
simParams sim = paramsSim->getSim();
+
for (int i = 0; i < (sim.years / sim.outIntOcc) + 1; i++) {
sum = ss = 0.0;
for (int rep = 0; rep < sim.reps; rep++) {
@@ -1195,15 +1218,20 @@ void Community::outOccSuit(bool view) {
if (sd > 0.0) sd = sqrt(sd);
else sd = 0.0;
se = sd / sqrt((double)(sim.reps));
+
outsuit << i * sim.outIntOcc << "\t" << mean << "\t" << se << endl;
- if (view) viewOccSuit(i * sim.outIntOcc, mean, se);
}
}
+// Close traits file
+bool Community::outTraitsFinishLandscape() {
+ return subComms[0]->outTraitsFinishLandscape();
+}
+
// Open traits file and write header record
-bool Community::outTraitsHeaders(Species* pSpecies, int landNr) {
- return subComms[0]->outTraitsHeaders(pLandscape, pSpecies, landNr);
+bool Community::outTraitsStartLandscape(Species* pSpecies, int landNr) {
+ return subComms[0]->outTraitsStartLandscape(pLandscape, pSpecies, landNr);
}
// Write records to traits file
@@ -1212,11 +1240,9 @@ only, this function relies on the fact that subcommunities are created in the sa
sequence as patches, which is in asecending order of x nested within descending
order of y
*/
-void Community::outTraits(traitCanvas tcanv, Species* pSpecies,
- int rep, int yr, int gen)
+void Community::outTraits(Species* pSpecies, int rep, int yr, int gen)
{
simParams sim = paramsSim->getSim();
- simView v = paramsSim->getViews();
landParams land = pLandscape->getLandParams();
traitsums* ts = 0;
traitsums sctraits;
@@ -1224,7 +1250,7 @@ void Community::outTraits(traitCanvas tcanv, Species* pSpecies,
// create array of traits means, etc., one for each row
ts = new traitsums[land.dimY];
for (int y = 0; y < land.dimY; y++) {
- for (int i = 0; i < NSEXES; i++) {
+ for (int i = 0; i < gMaxNbSexes; i++) {
ts[y].ninds[i] = 0;
ts[y].sumD0[i] = ts[y].ssqD0[i] = 0.0;
ts[y].sumAlpha[i] = ts[y].ssqAlpha[i] = 0.0;
@@ -1237,22 +1263,22 @@ void Community::outTraits(traitCanvas tcanv, Species* pSpecies,
ts[y].sumS0[i] = ts[y].ssqS0[i] = 0.0;
ts[y].sumAlphaS[i] = ts[y].ssqAlphaS[i] = 0.0;
ts[y].sumBetaS[i] = ts[y].ssqBetaS[i] = 0.0;
+ ts[y].sumGeneticFitness[i] = ts[y].ssqGeneticFitness[i] = 0.0;
}
}
}
- if (v.viewTraits
- || ((sim.outTraitsCells && yr >= sim.outStartTraitCell && yr % sim.outIntTraitCell == 0) ||
- (sim.outTraitsRows && yr >= sim.outStartTraitRow && yr % sim.outIntTraitRow == 0)))
+ if ((sim.outTraitsCells && yr >= sim.outStartTraitCell && yr % sim.outIntTraitCell == 0) ||
+ (sim.outTraitsRows && yr >= sim.outStartTraitRow && yr % sim.outIntTraitRow == 0))
{
// generate output for each sub-community (patch) in the community
int nsubcomms = (int)subComms.size();
for (int i = 1; i < nsubcomms; i++) { // // all except matrix sub-community
- sctraits = subComms[i]->outTraits(tcanv, pLandscape, rep, yr, gen, false);
+ sctraits = subComms[i]->outTraits(pLandscape, rep, yr, gen, false);
locn loc = subComms[i]->getLocn();
int y = loc.y;
if (sim.outTraitsRows && yr >= sim.outStartTraitRow && yr % sim.outIntTraitRow == 0)
{
- for (int s = 0; s < NSEXES; s++) {
+ for (int s = 0; s < gMaxNbSexes; s++) {
ts[y].ninds[s] += sctraits.ninds[s];
ts[y].sumD0[s] += sctraits.sumD0[s]; ts[y].ssqD0[s] += sctraits.ssqD0[s];
ts[y].sumAlpha[s] += sctraits.sumAlpha[s]; ts[y].ssqAlpha[s] += sctraits.ssqAlpha[s];
@@ -1265,6 +1291,7 @@ void Community::outTraits(traitCanvas tcanv, Species* pSpecies,
ts[y].sumS0[s] += sctraits.sumS0[s]; ts[y].ssqS0[s] += sctraits.ssqS0[s];
ts[y].sumAlphaS[s] += sctraits.sumAlphaS[s]; ts[y].ssqAlphaS[s] += sctraits.ssqAlphaS[s];
ts[y].sumBetaS[s] += sctraits.sumBetaS[s]; ts[y].ssqBetaS[s] += sctraits.ssqBetaS[s];
+ ts[y].sumGeneticFitness[s] += sctraits.sumGeneticFitness[s]; ts[y].ssqGeneticFitness[s] += sctraits.ssqGeneticFitness[s];
}
}
}
@@ -1284,8 +1311,8 @@ void Community::outTraits(traitCanvas tcanv, Species* pSpecies,
void Community::writeTraitsRows(Species* pSpecies, int rep, int yr, int gen, int y,
traitsums ts)
{
- emigRules emig = pSpecies->getEmig();
- trfrRules trfr = pSpecies->getTrfr();
+ emigRules emig = pSpecies->getEmigRules();
+ transferRules trfr = pSpecies->getTransferRules();
settleType sett = pSpecies->getSettle();
double mn, sd;
@@ -1347,7 +1374,7 @@ void Community::writeTraitsRows(Species* pSpecies, int rep, int yr, int gen, int
}
if (trfr.indVar) {
- if (trfr.moveModel) {
+ if (trfr.usesMovtProc) {
if (trfr.moveType == 2) { // CRW
// NB - CURRENTLY CANNOT BE SEX-DEPENDENT...
if (popsize > 0) mn = ts.sumStepL[0] / (double)popsize; else mn = 0.0;
@@ -1423,35 +1450,52 @@ void Community::writeTraitsRows(Species* pSpecies, int rep, int yr, int gen, int
if (popsize > 1) sd = ts.ssqBetaS[0] / (double)popsize - mn * mn; else sd = 0.0;
if (sd > 0.0) sd = sqrt(sd); else sd = 0.0;
outtraitsrows << "\t" << mn << "\t" << sd;
- // }
}
+ if (pSpecies->getNbGenLoadTraits() > 0) {
+ if (gMaxNbSexes > 1) {
+ if (ts.ninds[0] > 0) mn = ts.sumGeneticFitness[0] / (double)ts.ninds[0]; else mn = 0.0;
+ if (ts.ninds[0] > 1) sd = ts.ssqGeneticFitness[0] / (double)ts.ninds[0] - mn * mn; else sd = 0.0;
+ if (sd > 0.0) sd = sqrt(sd); else sd = 0.0;
+ outtraitsrows << "\t" << mn << "\t" << sd;
+ if (ts.ninds[1] > 0) mn = ts.sumGeneticFitness[1] / (double)ts.ninds[1]; else mn = 0.0;
+ if (ts.ninds[1] > 1) sd = ts.ssqGeneticFitness[1] / (double)ts.ninds[1] - mn * mn; else sd = 0.0;
+ if (sd > 0.0) sd = sqrt(sd); else sd = 0.0;
+ outtraitsrows << "\t" << mn << "\t" << sd;
+ }
+ else {
+ if (ts.ninds[0] > 0) mn = ts.sumGeneticFitness[0] / (double)ts.ninds[0]; else mn = 0.0;
+ if (ts.ninds[0] > 1) sd = ts.ssqGeneticFitness[0] / (double)ts.ninds[0] - mn * mn; else sd = 0.0;
+ if (sd > 0.0) sd = sqrt(sd); else sd = 0.0;
+ outtraitsrows << "\t" << mn << "\t" << sd;
+ }
+ }
outtraitsrows << endl;
}
-// Open trait rows file and write header record
-bool Community::outTraitsRowsHeaders(Species* pSpecies, int landNr) {
-
- if (landNr == -999) { // close file
- if (outtraitsrows.is_open()) outtraitsrows.close();
- outtraitsrows.clear();
- return true;
- }
+// Close trait rows file
+bool Community::outTraitsRowsFinishLandscape() {
+ if (outtraitsrows.is_open()) outtraitsrows.close();
+ outtraitsrows.clear();
+ return true;
+}
+// Open trait rows file and write header record
+bool Community::outTraitsRowsStartLandscape(Species* pSpecies, int landNr) {
string name;
- emigRules emig = pSpecies->getEmig();
- trfrRules trfr = pSpecies->getTrfr();
+ emigRules emig = pSpecies->getEmigRules();
+ transferRules trfr = pSpecies->getTransferRules();
settleType sett = pSpecies->getSettle();
simParams sim = paramsSim->getSim();
string DirOut = paramsSim->getDir(2);
if (sim.batchMode) {
name = DirOut
- + "Batch" + Int2Str(sim.batchNum) + "_"
- + "Sim" + Int2Str(sim.simulation) + "_Land" + Int2Str(landNr) + "_TraitsXrow.txt";
+ + "Batch" + to_string(sim.batchNum) + "_"
+ + "Sim" + to_string(sim.simulation) + "_Land" + to_string(landNr) + "_TraitsXrow.txt";
}
else {
- name = DirOut + "Sim" + Int2Str(sim.simulation) + "_TraitsXrow.txt";
+ name = DirOut + "Sim" + to_string(sim.simulation) + "_TraitsXrow.txt";
}
outtraitsrows.open(name.c_str());
@@ -1483,7 +1527,7 @@ bool Community::outTraitsRowsHeaders(Species* pSpecies, int landNr) {
}
}
if (trfr.indVar) {
- if (trfr.moveModel) {
+ if (trfr.usesMovtProc) {
if (trfr.moveType == 2) {
outtraitsrows << "\tmeanStepLength\tstdStepLength\tmeanRho\tstdRho";
}
@@ -1505,10 +1549,26 @@ bool Community::outTraitsRowsHeaders(Species* pSpecies, int landNr) {
}
if (sett.indVar) {
+ // if (sett.sexDep) {
+ // outtraitsrows << "\tF_meanS0\tF_stdS0\tM_meanS0\tM_stdS0";
+ // outtraitsrows << "\tF_meanAlphaS\tF_stdAlphaS\tM_meanAlphaS\tM_stdAlphaS";
+ // outtraitsrows << "\tF_meanBetaS\tF_stdBetaS\tM_meanBetaS\tM_stdBetaS";
+ // }
+ // else {
outtraitsrows << "\tmeanS0\tstdS0";
outtraitsrows << "\tmeanAlphaS\tstdAlphaS";
outtraitsrows << "\tmeanBetaS\tstdBetaS";
+ // }
}
+
+ if (pSpecies->getNbGenLoadTraits() > 0) {
+ if (gMaxNbSexes > 1) {
+ outtraitsrows << "\tF_meanProbViable\tF_stdProbViable\tM_meanProbViable\tM_stdProbViable";
+ }
+ else
+ outtraitsrows << "\tmeanProbViable\tstdProbViable";
+ }
+
outtraitsrows << endl;
return outtraitsrows.is_open();
@@ -1516,47 +1576,535 @@ bool Community::outTraitsRowsHeaders(Species* pSpecies, int landNr) {
}
#if RS_RCPP && !R_CMD
-Rcpp::IntegerMatrix Community::addYearToPopList(int rep, int yr) { // TODO: define new simparams to control start and interval of output
-
+Rcpp::IntegerMatrix Community::addYearToPopList(int rep, int yr, PopOutType type, int stage) { // TODO: define new simparams to control start and interval of output
+ /* Rcpp::Rcout << "Calling addYearToPopList: "
+ << "rep=" << rep
+ << " yr=" << yr
+ << " type=" << (int)type
+ << " stage=" << stage << endl;
+*/
landParams ppLand = pLandscape->getLandParams();
Rcpp::IntegerMatrix pop_map_year(ppLand.dimY, ppLand.dimX);
- intptr patch = 0;
Patch* pPatch = 0;
- intptr subcomm = 0;
- SubCommunity* pSubComm = 0;
+ SubCommunity* pSubComm = nullptr;
popStats pop;
pop.nInds = pop.nAdults = pop.nNonJuvs = 0;
for (int y = 0; y < ppLand.dimY; y++) {
for (int x = 0; x < ppLand.dimX; x++) {
- Cell* pCell = pLandscape->findCell(x, y);
+ Cell* pCell = pLandscape->findCell(x, y); //if (pLandscape->cells[y][x] == 0) {
if (pCell == 0) { // no-data cell
pop_map_year(ppLand.dimY - 1 - y, x) = NA_INTEGER;
}
else {
- patch = pCell->getPatch();
- if (patch == 0) { // matrix cell
+ pPatch = pCell->getPatch();
+ if (pPatch == nullptr) { // matrix cell
pop_map_year(ppLand.dimY - 1 - y, x) = 0;
}
else {
- pPatch = (Patch*)patch;
- subcomm = pPatch->getSubComm();
- if (subcomm == 0) { // check if sub-community exists
+ pSubComm = pPatch->getSubComm();
+ if (pSubComm == nullptr) { // check if sub-community exists
pop_map_year(ppLand.dimY - 1 - y, x) = 0;
}
else {
- pSubComm = (SubCommunity*)subcomm;
pop = pSubComm->getPopStats();
- pop_map_year(ppLand.dimY - 1 - y, x) = pop.nInds; // use indices like this because matrix gets transposed upon casting it into a raster on R-level
+
+ switch (type) {
+ case PopOutType::NInd:
+ pop_map_year(ppLand.dimY - 1 - y, x) = pop.nInds;
+ break;
+
+ case PopOutType::Stage:
+ pop_map_year(ppLand.dimY - 1 - y, x) = pSubComm->getNbInds(stage); // check if function is correct?
+ break;
+
+ case PopOutType::Juvs:
+ pop_map_year(ppLand.dimY - 1 - y, x) = pSubComm->getNbInds(0);
+ break;
+ }
+ // pop_map_year(ppLand.dimY - 1 - y, x) = pop.nInds; // use indices like this because matrix gets transposed upon casting it into a raster on R-level
+ //pop_map_year(ppLand.dimY-1-y,x) = pop.nAdults;
}
}
}
}
}
+ //list_outPop.push_back(pop_map_year, "rep" + std::to_string(rep) + "_year" + std::to_string(yr));
+ return pop_map_year;
+}
+
+// write a similar function for patch-based models;
+// Instead of a spatial x,y raster, the output should also be a Rcpp::IntegerMatrix with PatchID and population size (or stage-specific population size) for each patch.
+// The number of columns is determined by what the user specified in ReturnStages:
+// As default its 2 columns: 1st column is the patch ID, 2nd column is the total abundance.
+// Depending on the user specification the columns 3 to maximal (number of stages + 2)
+// can contain the abundance of each stage (e.g. column 3 is abundance of juveniles (stage 0), column 4 is abundance of stage 1 etc).
+// But only selected stages are included, so if the user only wants to output juveniles and adults,
+// then column 3 is abundance of juveniles (stage 0) and column 4 is abundance of adults (stage 1), and no other stages are included in the output.
+// After the runtime, the user can create a spatial raster in R by joining it with the patch coordinates.
+// be aware: the output is then not a spatial raster, but a table
+Rcpp::IntegerMatrix Community::addYearToPopListPatchBased(int rep, int yr, Rcpp::LogicalVector stages) {
+ /* Rcpp::Rcout << "Calling addYearToPopListPatchBased: "
+ << "rep=" << rep
+ << " yr=" << yr << endl;*/
+ int nrows=pLandscape->getPatchNbs().size();
+ int ncols = 2; // for patchID + total abundance
+ std::vector stageIndices;
+
+ for (int i = 0; i < stages.length(); i++) {
+ if (stages[i]) {
+ stageIndices.push_back(i);
+ }
+ }
+ ncols += stageIndices.size();
+
+
+ Rcpp::IntegerMatrix pop_map_year(nrows, ncols); // 2 columns: 1st column is the patch ID, 2nd column is the total abundance (or stage-specific abundance depending on user specification)
+ Patch* pPatch = nullptr;
+ SubCommunity* pSubComm = nullptr;
+ int currentRow = 0;
+
+ for (auto patchId : pLandscape->getPatchNbs()) {
+ // pPatch = pLandscape->findPatch(patchId);
+ // if (pPatch == nullptr) { // check if patch exists
+ // continue; // skip to next patch
+ // } else{
+ // pSubComm = pPatch->getSubComm();
+ // if (pSubComm == nullptr) { // check if sub-community exists
+ // pop = pSubComm->getPopStats();
+ // pop_map_year(patchId, 0) = patchId; // 1st column is patch ID
+ // pop_map_year(patchId, 1) = 0; // 2nd column is total abundance
+ // // additional columns for stage-specific abundances depending on user specification
+ // for(int i = 0; i < stages.length(); i++) {
+ // int ncol = 0;
+ // if(stages[i]) {
+ // pop_map_year(patchId, 2 + ncol) = 0; // all following columns are stage specific columns depending on the users specifications
+ // ncol++;
+ // }
+ // }
+ // } else {
+ // pop = pSubComm->getPopStats();
+ // pop_map_year(patchId, 0) = patchId; // 1st column is patch ID
+ // pop_map_year(patchId, 1) = pop.nInds; // 2nd column is total abundance
+ // // additional columns for stage-specific abundances depending on user specification
+ // for (int i = 0; i < stages.length(); i++) {
+ // int ncol = 0;
+ // if(stages[i]) {
+ // pop_map_year(patchId, 2 + ncol) = pSubComm->getNbInds(stages[i]); // all following columns are stage specific columns depending on the users specifications
+ // ncol++;
+ // }
+ // }
+ // }
+ // }
+ pop_map_year(currentRow, 0) = patchId; // 1st column: patch ID
+ // loop over ncols to fill default to 0
+ for(int i = 1; i < ncols; i++) {
+ pop_map_year(currentRow, i) = 0; // 2nd column: default total abundance
+ }
+
+ pPatch = pLandscape->findPatch(patchId);
+ if (pPatch != nullptr) { // Valid patch
+ pSubComm = pPatch->getSubComm();
+ if (pSubComm != nullptr) { // Valid sub-community
+ popStats pop = pSubComm->getPopStats();
+ pop_map_year(currentRow, 1) = pop.nInds; // Actual total abundance
+
+ for (int idx = 0; idx < stageIndices.size(); ++idx) {
+ int stage = stageIndices[idx];
+ pop_map_year(currentRow, 2 + idx) = pSubComm->getNbInds(stage);
+ }
+ }
+ }
+ currentRow++;
+ }
+
return pop_map_year;
+ // }
+ // return pop_map_year;
}
+
#endif
+bool Community::openOutGenesFile(const bool& isDiploid, const int landNr, const int rep)
+{
+ if (landNr == -999) { // close the file
+ if (ofsGenes.is_open()) {
+ ofsGenes.close();
+ ofsGenes.clear();
+ }
+ return true;
+ }
+ string name;
+ simParams sim = paramsSim->getSim();
+
+ if (sim.batchMode) {
+ name = paramsSim->getDir(2)
+ + "Batch" + to_string(sim.batchNum) + "_"
+ + "Sim" + to_string(sim.simulation) + "_Land"
+ + to_string(landNr) + "_Rep" + to_string(rep)
+ + "_geneValues.txt";
+ }
+ else {
+ name = paramsSim->getDir(2)
+ + "Sim" + to_string(sim.simulation) + "_Land"
+ + to_string(landNr) + "_Rep" + to_string(rep)
+ + "_geneValues.txt";
+ }
+
+ ofsGenes.open(name.c_str());
+ ofsGenes << "Year\tGeneration\tIndID\ttraitType\tlocusPosition"
+ << "\talleleValueA\tdomCoefA";
+ if (isDiploid) ofsGenes << "\talleleValueB\tdomCoefB";
+ ofsGenes << endl;
+
+ return ofsGenes.is_open();
+}
+
+void Community::outputGeneValues(const int& year, const int& gen, Species* pSpecies) {
+ if (!ofsGenes.is_open())
+ throw runtime_error("Could not open output gene values file.");
+
+ const set patchList = pSpecies->getSamplePatches();
+ for (int patchId : patchList) {
+ const auto patch = pLandscape->findPatch(patchId);
+ if (patch == 0) {
+ throw runtime_error("Sampled patch does not exist");
+ }
+ const auto pPop = patch->getPopn(pSpecies);
+ if (pPop != 0) {
+ pPop->outputGeneValues(ofsGenes, year, gen);
+ }
+ }
+}
+
+// ----------------------------------------------------------------------------------------
+// Sample individuals from sample patches
+// ----------------------------------------------------------------------------------------
+
+void Community::sampleIndividuals(Species* pSpecies) {
+
+ const set patchList = pSpecies->getSamplePatches();
+ string nbIndsToSample = pSpecies->getNIndsToSample();
+ const set stagesToSampleFrom = pSpecies->getStagesToSample();
+
+ for (int patchId : patchList) {
+ const auto patch = pLandscape->findPatch(patchId);
+ if (patch == 0) {
+ throw runtime_error("Can't sample individuals: patch" + to_string(patchId) + "doesn't exist.");
+ }
+ auto pPop = patch->getPopn(pSpecies);
+ if (pPop != 0) {
+ pPop->sampleIndsWithoutReplacement(nbIndsToSample, stagesToSampleFrom);
+ }
+ }
+}
+
+// ----------------------------------------------------------------------------------------
+// Open population level Fstat output file
+// ----------------------------------------------------------------------------------------
+
+bool Community::openNeutralOutputFile(Species* pSpecies, int landNr)
+{
+ if (landNr == -999) { // close the file
+ if (outwcfstat.is_open()) outwcfstat.close();
+ outwcfstat.clear();
+ return true;
+ }
+
+ string name;
+ simParams sim = paramsSim->getSim();
+
+ if (sim.batchMode) {
+ name = paramsSim->getDir(2)
+ + "Batch" + to_string(sim.batchNum) + "_"
+ + "Sim" + to_string(sim.simulation) + "_Land"
+ + to_string(landNr)
+ + "_neutralGenetics.txt";
+ }
+ else {
+ name = paramsSim->getDir(2) + "Sim" + to_string(sim.simulation) + "_neutralGenetics.txt";
+ }
+ outwcfstat.open(name.c_str());
+ outwcfstat << "Rep\tYear\tRepSeason\tnExtantPatches\tnIndividuals\tFstWC\tFisWC\tFitWC\tmeanAllelePerLocus\tmeanAllelePerLocusPatches\tmeanFixedLoci\tmeanFixedLociPatches\tmeanObHeterozygosity";
+ outwcfstat << endl;
+
+ return outwcfstat.is_open();
+}
+
+// ----------------------------------------------------------------------------------------
+// open per locus WC fstat using MS approach, this will output MS calculated FIS, FIT, FST
+// in general population neutral genetics output file
+// ----------------------------------------------------------------------------------------
+
+bool Community::openPerLocusFstFile(Species* pSpecies, Landscape* pLandscape, const int landNr, const int rep)
+{
+ set patchList;
+ string samplingOpt = paramsSim->getSim().patchSamplingOption;
+ if (samplingOpt == "list"
+ || samplingOpt == "random"
+ || (samplingOpt == "all" && !pLandscape->getLandParams().dynamic)
+ ) // list of patches always the same
+ patchList = pSpecies->getSamplePatches();
+ else { // random_occupied or all with dynamic landscape
+ // then sampled patches may change every year,
+ // so produce an entry for all patches
+ patchList = pLandscape->getPatchNbs();
+ }
+
+ if (landNr == -999) { // close the file
+ if (outperlocusfstat.is_open()) outperlocusfstat.close();
+ outperlocusfstat.clear();
+ return true;
+ }
+
+ string name;
+ simParams sim = paramsSim->getSim();
+
+ if (sim.batchMode) {
+ name = paramsSim->getDir(2)
+ + "Batch" + to_string(sim.batchNum) + "_"
+ + "Sim" + to_string(sim.simulation) + "_Land"
+ + to_string(landNr) + "_Rep"
+ + to_string(rep)
+ + "_perLocusNeutralGenetics.txt";
+ }
+ else {
+ name = paramsSim->getDir(2) + "Sim" + to_string(sim.simulation) + "_Rep" + to_string(rep) + "_perLocusNeutralGenetics.txt";
+ }
+
+ outperlocusfstat.open(name.c_str());
+ outperlocusfstat << "Year\tRepSeason\tLocus\tFst\tFis\tFit\tHet";
+ for (int patchId : patchList) {
+ outperlocusfstat << "\tpatch_" + to_string(patchId) + "_Het";
+ }
+ outperlocusfstat << endl;
+
+ return outperlocusfstat.is_open();
+}
+
+// ----------------------------------------------------------------------------------------
+// open pairwise fst file
+// ----------------------------------------------------------------------------------------
+
+bool Community::openPairwiseFstFile(Species* pSpecies, Landscape* pLandscape, const int landNr, const int rep) {
+
+ const set patchList = pSpecies->getSamplePatches();
+
+ if (landNr == -999) { // close the file
+ if (outpairwisefst.is_open()) outpairwisefst.close();
+ outpairwisefst.clear();
+ return true;
+ }
+
+ string name;
+ simParams sim = paramsSim->getSim();
+
+ if (sim.batchMode) {
+ name = paramsSim->getDir(2)
+ + "Batch" + to_string(sim.batchNum) + "_"
+ + "Sim" + to_string(sim.simulation) + "_Land"
+ + to_string(landNr) + "_Rep"
+ + to_string(rep)
+ + "_pairwisePatchNeutralGenetics.txt";
+ }
+ else {
+ name = paramsSim->getDir(2) + "Sim" + to_string(sim.simulation) + "_Rep" + to_string(rep) + "_pairwisePatchNeutralGenetics.txt";
+ }
+ outpairwisefst.open(name.c_str());
+ outpairwisefst << "Year\tRepSeason\tpatchA\tpatchA_x\tpatchA_y\tpatchB\tpatchB_x\tpatchB_y\tFst";
+ outpairwisefst << endl;
+
+ return outpairwisefst.is_open();
+}
+
+// ----------------------------------------------------------------------------------------
+// Write population level FST results file
+// ----------------------------------------------------------------------------------------
+
+void Community::writeNeutralOutputFile(int rep, int yr, int gen) {
+
+ outwcfstat << rep << "\t" << yr << "\t" << gen << "\t";
+ outwcfstat << pNeutralStatistics->getNbPopulatedSampledPatches()
+ << "\t" << pNeutralStatistics->getTotalNbSampledInds() << "\t";
+
+
+ outwcfstat << pNeutralStatistics->getFstWC() << "\t"
+ << pNeutralStatistics->getFisWC() << "\t"
+ << pNeutralStatistics->getFitWC() << "\t";
+
+
+
+ outwcfstat << pNeutralStatistics->getMeanNbAllPerLocus() << "\t"
+ << pNeutralStatistics->getMeanNbAllPerLocusPerPatch() << "\t"
+ << pNeutralStatistics->getTotalFixdAlleles() << "\t"
+ << pNeutralStatistics->getMeanFixdAllelesPerPatch() << "\t"
+ << pNeutralStatistics->getHo();
+
+ outwcfstat << endl;
+}
+
+// ----------------------------------------------------------------------------------------
+// Write per locus FST results file
+// ----------------------------------------------------------------------------------------
+
+void Community::writePerLocusFstatFile(Species* pSpecies, const int yr, const int gen, const int nLoci, set const& patchList)
+{
+ string samplingOpt = paramsSim->getSim().patchSamplingOption;
+ bool samplingFixed = samplingOpt == "list"
+ || samplingOpt == "random"
+ || (samplingOpt == "all" && !pLandscape->getLandParams().dynamic);
+
+ const set positions = pSpecies->getSpTrait(NEUTRAL)->getGenePositions();
+
+ int thisLocus = 0;
+ for (int position : positions) {
+
+ outperlocusfstat << yr << "\t"
+ << gen << "\t"
+ << position << "\t";
+ outperlocusfstat << pNeutralStatistics->getPerLocusFst(thisLocus) << "\t"
+ << pNeutralStatistics->getPerLocusFis(thisLocus) << "\t"
+ << pNeutralStatistics->getPerLocusFit(thisLocus) << "\t"
+ << pNeutralStatistics->getPerLocusHo(thisLocus);
+
+ if (samplingFixed) { // then safe to output sampled patches in order
+ for (int patchId : patchList) {
+ float het = getPatchHet(pSpecies, patchId, thisLocus);
+ if (het < 0) // patch empty
+ outperlocusfstat << "\t" << "NA";
+ else outperlocusfstat << "\t" << het;
+ }
+ }
+ else { // sampling may change between generations, must produce output for all patches in Landscape
+ for (auto patchId : pLandscape->getPatchNbs()) {
+ if (patchList.find(patchId) != patchList.end()) {
+ float het = getPatchHet(pSpecies, patchId, thisLocus);
+ if (het < 0) // patch empty
+ outperlocusfstat << "\t" << "NA";
+ else outperlocusfstat << "\t" << het;
+ }
+ else { // patch not sampled
+ outperlocusfstat << "\t" << "NA";
+ }
+ }
+ }
+ ++thisLocus;
+ outperlocusfstat << endl;
+ }
+}
+
+// Calculate the observed heterozygosity (Ho) for a patch
+// = number of heterozygous individuals at this locus / nb individuals in patch
+float Community::getPatchHet(Species* pSpecies, int patchId, int whichLocus) const {
+ const auto patch = pLandscape->findPatch(patchId);
+ const auto pPop = patch->getPopn(pSpecies);
+ int nAlleles = pSpecies->getSpTrait(NEUTRAL)->getNbNeutralAlleles();
+ int popSize = 0;
+ float het = 0;
+ if (pPop != 0) {
+ popSize = pPop->sampleSize();
+ if (popSize == 0) return -1.0;
+ else {
+ for (int a = 0; a < nAlleles; ++a) {
+ het += static_cast(pPop->getHeteroTally(whichLocus, a));
+ }
+ het /= popSize;
+ return het;
+ }
+ }
+ else return -1.0;
+}
+
+
+// ----------------------------------------------------------------------------------------
+// Write pairwise FST results file
+// ----------------------------------------------------------------------------------------
+void Community::writePairwiseFstFile(Species* pSpecies, const int yr, const int gen, set const& patchList) {
+
+ const int nPatches = static_cast(patchList.size());
+ // Convert set to vector for index-based access
+ vector patchVect;
+ copy(patchList.begin(), patchList.end(), back_inserter(patchVect));
+
+ for (int i = 0; i < nPatches; ++i)
+ {
+ const auto patchA = pLandscape->findPatch(patchVect[i]);
+
+ for (int j = i; j < nPatches; ++j)
+ {
+
+ const auto patchB = pLandscape->findPatch(patchVect[j]);
+
+ outpairwisefst << yr << "\t"
+ << gen << "\t"
+ << patchVect[i] << "\t"
+ << patchA->getSubComm()->getLocn().x << "\t"
+ << patchA->getSubComm()->getLocn().y << "\t"
+ << patchVect[j] << "\t"
+ << patchB->getSubComm()->getLocn().x << "\t"
+ << patchB->getSubComm()->getLocn().y << "\t"
+ << pNeutralStatistics->getPairwiseFst(i, j)
+ << "\n";
+ }
+ }
+
+
+}
+
+
+// ----------------------------------------------------------------------------------------
+// Output and calculate neutral statistics
+// ----------------------------------------------------------------------------------------
+
+
+void Community::calculateNeutralGenetics(Species* pSpecies, int rep, int yr, int gen, bool outPairwiseFst, int outputPairwiseFstStart, int outputPairwiseFstInterval,
+ bool outputGlobalFst, int outputGlobalFstStart, int outputGlobalFstInterval, bool outputPerLocusFst) {
+
+ const int maxNbNeutralAlleles = pSpecies->getSpTrait(NEUTRAL)->getNbNeutralAlleles();
+ const int nLoci = (int)pSpecies->getNPositionsForTrait(NEUTRAL);
+ const set patchList = pSpecies->getSamplePatches();
+ int nInds = 0, nbPops = 0;
+
+ for (int patchId : patchList) {
+ const auto patch = pLandscape->findPatch(patchId);
+ if (patch == 0) {
+ throw runtime_error("Sampled patch does not exist");
+ }
+ const auto pPop = patch->getPopn(pSpecies);
+ if (pPop != 0) { // empty patches do not contribute
+ nInds += pPop->sampleSize();
+ nbPops++;
+ }
+ }
+
+ if (pNeutralStatistics == 0)
+ pNeutralStatistics = make_unique(patchList.size(), nLoci);
+
+ pNeutralStatistics->updateAllNeutralTables(pSpecies, pLandscape, patchList);
+ pNeutralStatistics->calculateHo(patchList, nInds, nLoci, pSpecies, pLandscape);
+ pNeutralStatistics->calculatePerLocusHo(patchList, nInds, nLoci, pSpecies, pLandscape);
+ pNeutralStatistics->calcAllelicDiversityMetrics(patchList, nInds, pSpecies, pLandscape);
+
+ if (outPairwiseFst) {
+ pNeutralStatistics->calculatePairwiseFst(patchList, nLoci, maxNbNeutralAlleles, pSpecies, pLandscape);
+
+ if (yr >= outputPairwiseFstStart && yr % outputPairwiseFstInterval == 0) {
+ writePairwiseFstFile(pSpecies, yr, gen, patchList);
+ }
+ }
+ if (outputGlobalFst) {
+ pNeutralStatistics->calculateFstatWC(patchList, nInds, nLoci, maxNbNeutralAlleles, pSpecies, pLandscape, false);
+
+ if (yr >= outputGlobalFstStart && yr % outputGlobalFstInterval == 0) {
+ writeNeutralOutputFile(rep, yr, gen);
+ if (outputPerLocusFst)
+ writePerLocusFstatFile(pSpecies, yr, gen, nLoci, patchList);
+ }
+ }
+
+}
+
+
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
diff --git a/Community.h b/Community.h
index c55b333..f102fd1 100644
--- a/Community.h
+++ b/Community.h
@@ -1,25 +1,25 @@
/*----------------------------------------------------------------------------
- *
- * Copyright (C) 2020 Greta Bocedi, Stephen C.F. Palmer, Justin M.J. Travis, Anne-Kathleen Malchow, Damaris Zurell
- *
+ *
+ * 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 .
- *
+ *
--------------------------------------------------------------------------*/
-
-
+
+
/*------------------------------------------------------------------------------
RangeShifter v2.0 Community
@@ -35,9 +35,9 @@ Optionally, the Community maintains a record of the occupancy of suitable cells
or patches during the course of simulation of multiple replicates.
For full details of RangeShifter, please see:
-Bocedi G., Palmer S.C.F., Peer G., Heikkinen R.K., Matsinos Y.G., Watts K.
+ Bocedi G., Palmer S.C.F., Pe’er G., Heikkinen R.K., Matsinos Y.G., Watts K.
and Travis J.M.J. (2014). RangeShifter: a platform for modelling spatial
-eco-evolutionary dynamics and species responses to environmental changes.
+ eco-evolutionary dynamics and species’ responses to environmental changes.
Methods in Ecology and Evolution, 5, 388-396. doi: 10.1111/2041-210X.12162
Authors: Greta Bocedi & Steve Palmer, University of Aberdeen
@@ -51,6 +51,8 @@ Last updated: 25 June 2021 by Anne-Kathleen Malchow
#include
#include
+#include
+#include
using namespace std;
#include "SubCommunity.h"
@@ -58,6 +60,7 @@ using namespace std;
#include "Patch.h"
#include "Cell.h"
#include "Species.h"
+#include "NeutralStatsManager.h"
//---------------------------------------------------------------------------
struct commStats {
@@ -65,6 +68,14 @@ int ninds,nnonjuvs,suitable,occupied;
int minX,maxX,minY,maxY;
};
+#if RS_RCPP// For raster output only: which type of population output should be stored?
+enum class PopOutType {
+ NInd, // total abundance
+ Stage, // specific stages
+ Juvs // juvenile stage
+};
+#endif
+
class Community {
public:
@@ -95,15 +106,14 @@ class Community {
);
#endif // SEASONAL || RS_RCPP
- void survival(
- short, // part: 0 = determine survival & development,
- // 1 = apply survival changes to the population
+ void survival0( // Determine survival & development
short, // option0: 0 = stage 0 (juveniles) only )
// 1 = all stages ) used by part 0 only
// 2 = stage 1 and above (all non-juvs) )
short // option1: 0 - development only (when survival is annual)
// 1 - development and survival
);
+ void survival1(); // Apply survival changes to the population
void ageIncrement(void);
int totalInds(void);
Population* findPop( // Find the population of a given species in a given patch
@@ -123,9 +133,10 @@ class Community {
int // no. of rows (as above)
);
- bool outRangeHeaders( // Open range file and write header record
+ bool outRangeFinishLandscape(); // Close range file
+ bool outRangeStartLandscape( // Open range file and write header record
Species*, // pointer to Species
- int // Landscape number (-999 to close the file)
+ int // Landscape number
);
void outRange( // Write record to range file
Species*, // pointer to Species
@@ -133,9 +144,9 @@ class Community {
int, // year
int // generation
);
- bool outPopHeaders( // Open population file and write header record
- Species*, // pointer to Species
- int // option: -999 to close the file
+ bool outPopFinishLandscape(); // Close population file
+ bool outPopStartLandscape( // Open population file and write header record
+ Species* // pointer to Species
);
void outPop( // Write records to population file
int, // replicate
@@ -143,46 +154,33 @@ class Community {
int // generation
);
- void outInds( // Write records to individuals file
+ void outIndsFinishReplicate(); // Close individuals file
+ void outIndsStartReplicate( // Open individuals file and write header record
int, // replicate
- int, // year
- int, // generation
- int // Landscape number (>= 0 to open the file, -999 to close the file
- // -1 to write data records)
+ int // Landscape number
);
- void outGenetics( // Write records to genetics file
+ void outIndividuals( // Write records to individuals file
int, // replicate
int, // year
- int, // generation
- int // Landscape number (>= 0 to open the file, -999 to close the file
- // -1 to write data records)
+ int // generation
);
+ // Close occupancy file
+ bool outOccupancyFinishLandscape();
// Open occupancy file, write header record and set up occupancy array
- bool outOccupancyHeaders(
- int // option: -999 to close the file
- );
+ bool outOccupancyStartLandscape();
void outOccupancy(void);
- void outOccSuit(
- bool // TRUE if occupancy graph is to be viewed on screen
- );
- void viewOccSuit( // Update the occupancy graph on the screen
- // NULL for the batch version
- int, // year
- double, // mean occupancy
- double // standard error of occupancy
- );
- bool outTraitsHeaders( // Open traits file and write header record
+ void outOccSuit();
+ bool outTraitsFinishLandscape(); // Close traits file
+ bool outTraitsStartLandscape( // Open traits file and write header record
Species*, // pointer to Species
- int // Landscape number (-999 to close the file)
+ int // Landscape number
);
- bool outTraitsRowsHeaders( // Open trait rows file and write header record
+ bool outTraitsRowsFinishLandscape(); // Close trait rows file
+ bool outTraitsRowsStartLandscape( // Open trait rows file and write header record
Species*, // pointer to Species
- int // Landscape number (-999 to close the file)
+ int // Landscape number
);
void outTraits( // Write records to traits file
- traitCanvas,// pointers to canvases for drawing variable traits
- // see SubCommunity.h
- // in the batch version, these are replaced by integers set to zero
Species*, // pointer to Species
int, // replicate
int, // year
@@ -196,23 +194,42 @@ class Community {
int, // row number (Y cell co-ordinate)
traitsums // structure holding sums of trait genes for dispersal (see Population.h)
);
- void draw( // Draw the Community on the landscape map and optionally save the map
- // NULL for the batch version
- int, // replicate
- int, // year
- int, // generation
- int // Landscape number
- );
#if RS_RCPP && !R_CMD
- Rcpp::IntegerMatrix addYearToPopList(int,int);
+ Rcpp::IntegerMatrix addYearToPopList(int,int,PopOutType,int);
+
+ Rcpp::IntegerMatrix addYearToPopListPatchBased(int,int,Rcpp::LogicalVector);
#endif
+ //sample individuals for genetics (or could be used for anything)
+ void sampleIndividuals(Species* pSpecies);
+
+ bool openOutGenesFile(const bool& isDiploid, const int landNr, const int rep);
+ void outputGeneValues(const int& year, const int& gen, Species* pSpecies);
+
+ //control neutral stat output
+
+ void calculateNeutralGenetics(Species* pSpecies, int rep, int yr, int gen, bool outPairwiseFst, int outputPairwiseFstStart, int outputPairwiseFstInterval,
+ bool outputGlobalFst, int outputGlobalFstStart, int outputGlobalFstInterval, bool outputPerLocusFst);
+
+
+ //file openers
+ bool openNeutralOutputFile(Species* pSpecies, const int landNr);
+ bool openPerLocusFstFile(Species* pSpecies, Landscape* pLandscape, const int landNr, const int rep);
+ bool openPairwiseFstFile(Species* pSpecies, Landscape* pLandscape, const int landNr, const int rep);
+
+ //file writers
+ void writeNeutralOutputFile(int rep, int yr, int gen);
+ void writePerLocusFstatFile(Species* pSpecies, const int yr, const int gen, const int nLoci, set const& patchList);
+ void writePairwiseFstFile(Species* pSpecies, const int yr, const int gen, set const& patchList);
+ float getPatchHet(Species* pSpecies, int patchId, int whichLocus) const;
private:
Landscape *pLandscape;
int indIx; // index used to apply initial individuals
float **occSuit; // occupancy of suitable cells / patches
std::vector subComms;
+ //below won't work for multispecies
+ unique_ptr pNeutralStatistics;
};
extern paramSim *paramsSim;
diff --git a/DispersalTrait.cpp b/DispersalTrait.cpp
new file mode 100644
index 0000000..f0a805d
--- /dev/null
+++ b/DispersalTrait.cpp
@@ -0,0 +1,489 @@
+/*----------------------------------------------------------------------------
+ *
+ * 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 .
+ *
+ * File Created by Roslyn Henry March 2023. Code adapted from NEMO (https://nemo2.sourceforge.io/)
+ --------------------------------------------------------------------------*/
+
+#include "DispersalTrait.h"
+
+// ----------------------------------------------------------------------------------------
+// Initialisation constructor
+// Called when initialising community
+// Sets up initial values, and immutable attributes (distributions and parameters)
+// that are defined at the species-level
+// ----------------------------------------------------------------------------------------
+DispersalTrait::DispersalTrait(SpeciesTrait* P)
+{
+ pSpeciesTrait = P;
+ ExpressionType expressionType = pSpeciesTrait->getExpressionType();
+
+ if (!pSpeciesTrait->isInherited()) // there is a trait for individual variation but this isn't inherited variation it's sampled from initial distribution
+ _inherit_func_ptr = &DispersalTrait::reInitialiseGenes;
+ else {
+ _inherit_func_ptr = (pSpeciesTrait->getPloidy() == 1) ? &DispersalTrait::inheritHaploid : &DispersalTrait::inheritDiploid; //this could be changed if we wanted some alternative form of inheritance
+
+ DistributionType mutationDistribution = pSpeciesTrait->getMutationDistribution();
+ map mutationParameters = pSpeciesTrait->getMutationParameters();
+
+ // Set mutation parameters
+ switch (mutationDistribution) {
+ case UNIFORM:
+ {
+ if (mutationParameters.count(MAX) != 1)
+ throw logic_error("mutation uniform distribution parameter must contain max value (e.g. max= ) \n");
+ if (mutationParameters.count(MIN) != 1)
+ throw logic_error("mutation uniform distribution parameter must contain min value (e.g. min= ) \n");
+ _mutate_func_ptr = &DispersalTrait::mutateUniform;
+ break;
+ }
+ case NORMAL:
+ {
+ if (mutationParameters.count(MEAN) != 1)
+ throw logic_error("mutation distribution set to normal so parameters must contain mean value (e.g. mean= ) \n");
+ if (mutationParameters.count(SD) != 1)
+ throw logic_error("mutation distribution set to normal so parameters must contain sdev value (e.g. sdev= ) \n");
+ _mutate_func_ptr = &DispersalTrait::mutateNormal;
+ break;
+ }
+ default:
+ {
+ throw logic_error("wrong parameter value for mutation model, must be uniform/normal \n"); //unless want to add gamma or negative exp
+ break;
+ }
+ }
+ }
+
+ // Set initialisation parameters
+ DistributionType initialDistribution = pSpeciesTrait->getInitialDistribution();
+ map initialParameters = pSpeciesTrait->getInitialParameters();
+ switch (initialDistribution) {
+ case UNIFORM:
+ {
+ if (initialParameters.count(MAX) != 1)
+ throw logic_error("initial uniform distribution parameter must contain max value (e.g. max= ) \n");
+ if (initialParameters.count(MIN) != 1)
+ throw logic_error("initial uniform distribution parameter must contain min value (e.g. min= ) \n");
+ float maxD = initialParameters.find(MAX)->second;
+ float minD = initialParameters.find(MIN)->second;
+ initialiseUniform(minD, maxD);
+ break;
+ }
+ case NORMAL:
+ {
+ if (initialParameters.count(MEAN) != 1)
+ throw logic_error("initial normal distribution parameter must contain mean value (e.g. mean= ) \n");
+ if (initialParameters.count(SD) != 1)
+ throw logic_error("initial normal distribution parameter must contain sdev value (e.g. sdev= ) \n");
+ float mean = initialParameters.find(MEAN)->second;
+ float sd = initialParameters.find(SD)->second;
+ initialiseNormal(mean, sd);
+ break;
+ }
+ default:
+ {
+ throw logic_error("wrong parameter value for parameter \"initialisation of dispersal traits\", must be uniform/normal \n");
+ break;
+ }
+ }
+
+ // Set expression mode parameters
+ switch (expressionType) {
+ case AVERAGE:
+ {
+ _express_func_ptr = &DispersalTrait::expressAverage;
+ break;
+ }
+ case ADDITIVE:
+ {
+ _express_func_ptr = &DispersalTrait::expressAdditive;
+ break;
+ }
+ default:
+ {
+ throw logic_error("wrong parameter value for parameter \"expression of dispersal trait\", must be average/additive \n");
+ break;
+ }
+ }
+}
+
+// ----------------------------------------------------------------------------------------
+// Inheritance constructor
+// Copies immutable features from a parent trait
+// Only called via clone()
+// ----------------------------------------------------------------------------------------
+DispersalTrait::DispersalTrait(const DispersalTrait& T) :
+ pSpeciesTrait(T.pSpeciesTrait),
+ _mutate_func_ptr(T._mutate_func_ptr),
+ _inherit_func_ptr(T._inherit_func_ptr),
+ _express_func_ptr(T._express_func_ptr) {}
+
+// ----------------------------------------------------------------------------------------
+// Sample and apply mutations from a uniform distribution
+//
+// Mutations drawn only for existing positions,
+// that is no new genes are created during simulation
+// ----------------------------------------------------------------------------------------
+void DispersalTrait::mutateUniform()
+{
+ const int positionsSize = pSpeciesTrait->getPositionsSize();
+ const auto& genePositions = pSpeciesTrait->getGenePositions();
+ const short ploidy = pSpeciesTrait->getPloidy();
+ const float mutationRate = pSpeciesTrait->getMutationRate();
+ float newAlleleVal;
+
+ auto rng = pRandom->getRNG();
+
+ map mutationParameters = pSpeciesTrait->getMutationParameters();
+ float maxD = mutationParameters.find(MAX)->second;
+ float minD = mutationParameters.find(MIN)->second;
+
+ for (int p = 0; p < ploidy; p++) {
+
+ unsigned int NbMut = pRandom->Binomial(positionsSize, mutationRate);
+
+ if (NbMut > 0) {
+ vector mutationPositions;
+ sample(genePositions.begin(), genePositions.end(), std::back_inserter(mutationPositions),
+ NbMut, rng);
+
+ for (int m : mutationPositions) {
+ auto it = genes.find(m);
+ if (it == genes.end())
+ throw runtime_error("Locus sampled for mutation doesn't exist.");
+ float currentAlleleVal = it->second[p].get()->getAlleleValue();//current
+ newAlleleVal = pRandom->FRandom(minD, maxD) + currentAlleleVal;
+ it->second[p] = make_shared(newAlleleVal, dispDominanceFactor);
+ }
+ }
+ }
+}
+
+// ----------------------------------------------------------------------------------------
+// Sample and apply mutations from a normal distribution
+// Mutations drawn only for existing positions,
+// that is no new genes are created during simulation
+// ----------------------------------------------------------------------------------------
+void DispersalTrait::mutateNormal()
+{
+ const int positionsSize = pSpeciesTrait->getPositionsSize();
+ const auto& genePositions = pSpeciesTrait->getGenePositions();
+ const short ploidy = pSpeciesTrait->getPloidy();
+ const float mutationRate = pSpeciesTrait->getMutationRate();
+
+ auto rng = pRandom->getRNG();
+
+ const map mutationParameters = pSpeciesTrait->getMutationParameters();
+ const float mean = mutationParameters.find(MEAN)->second;
+ const float sd = mutationParameters.find(SD)->second;
+ float newAlleleVal;
+
+ for (int p = 0; p < ploidy; p++) {
+
+ unsigned int NbMut = pRandom->Binomial(positionsSize, mutationRate);
+
+ if (NbMut > 0) {
+ vector mutationPositions;
+ sample(genePositions.begin(), genePositions.end(), std::back_inserter(mutationPositions),
+ NbMut, rng);
+
+ for (int m : mutationPositions) {
+ auto it = genes.find(m);
+ if (it == genes.end())
+ throw runtime_error("Locus sampled for mutation doesn't exist.");
+ float currentAlleleVal = it->second[p].get()->getAlleleValue(); //current
+ newAlleleVal = pRandom->Normal(mean, sd) + currentAlleleVal;
+ it->second[p] = make_shared(newAlleleVal, dispDominanceFactor);
+ }
+ }
+ }
+}
+
+// ----------------------------------------------------------------------------------------
+// Wrapper to inheritance function
+// ----------------------------------------------------------------------------------------
+void DispersalTrait::inheritGenes(const bool& fromMother, QuantitativeTrait* parentTrait, set const& recomPositions, int startingChromosome)
+{
+ auto parentCast = dynamic_cast(parentTrait); // must convert QuantitativeTrait to DispersalTrait
+ const auto& parent_seq = parentCast->getGenes();
+ (this->*_inherit_func_ptr)(fromMother, parent_seq, recomPositions, startingChromosome);
+}
+
+// ----------------------------------------------------------------------------------------
+// Inheritance for diploid, sexual species
+// Called once for each parent.
+// Pass the correct parental strand, resolving crossing-overs
+// after each recombinant site e.g. if parent genotype is
+// 0000
+// 1111
+// and position 2 is selected to recombine, then offspring inherits
+// 0001
+// Assumes mother genes are inherited first.
+// ----------------------------------------------------------------------------------------
+void DispersalTrait::inheritDiploid(const bool& fromMother, map>> const& parentGenes, set const& recomPositions, int parentChromosome) {
+
+ const int lastPosition = parentGenes.rbegin()->first;
+ auto recomIt = recomPositions.lower_bound(parentGenes.begin()->first);
+ // If no recombination sites, only breakpoint is last position
+ // i.e., no recombination occurs
+ int nextBreakpoint = recomIt == recomPositions.end() ? lastPosition : *recomIt;
+
+ // Is the first parent gene position already recombinant?
+ auto distance = std::distance(recomPositions.begin(), recomIt);
+ if (distance % 2 != 0) // odd positions = switch, even = switch back
+ parentChromosome = 1 - parentChromosome; //switch chromosome
+
+ for (auto const& [locus, allelePair] : parentGenes) {
+
+ // Switch chromosome if locus is past recombination site
+ while (locus > nextBreakpoint) {
+ parentChromosome = 1 - parentChromosome;
+ std::advance(recomIt, 1); // go to next recombination site
+ nextBreakpoint = recomIt == recomPositions.end() ? lastPosition : *recomIt;
+ }
+
+ if (locus <= nextBreakpoint) {
+ auto& parentAllele = allelePair[parentChromosome];
+ auto itGene = genes.find(locus);
+ if (itGene == genes.end()) {
+ // locus does not exist yet, create and initialise it
+ if (!fromMother) throw runtime_error("Father-inherited locus does not exist.");
+ vector> newAllelePair(2);
+ newAllelePair[sex_t::FEM] = parentAllele;
+ genes.insert(make_pair(locus, newAllelePair));
+ }
+ else { // father, locus already exists
+ if (fromMother) throw runtime_error("Mother-inherited locus already exists.");
+ itGene->second[sex_t::MAL] = parentAllele;
+ }
+ }
+ }
+}
+
+// ----------------------------------------------------------------------------------------
+// Inheritance for haploid, asexual species
+// Simply pass down parent genes
+// Arguments are still needed to match overloaded function in base class
+// ----------------------------------------------------------------------------------------
+void DispersalTrait::inheritHaploid(const bool& fromMother, map>> const& parentGenes, set const& recomPositions, int parentChromosome)
+{
+ genes = parentGenes;
+}
+
+// ----------------------------------------------------------------------------------------
+// Non-inheritance
+// For cases where isInherited option is turned off
+// In this case, offspring alleles are populated using the initialise functions
+// Arguments are still needed to match overloaded function in base class
+// ----------------------------------------------------------------------------------------
+void DispersalTrait::reInitialiseGenes(const bool& fromMother, map>> const& parentGenes, set const& recomPositions, int parentChromosome)
+{
+ DistributionType initialDistribution = pSpeciesTrait->getInitialDistribution();
+ map initialParameters = pSpeciesTrait->getInitialParameters();
+
+ switch (initialDistribution) {
+ case UNIFORM:
+ {
+ if (initialParameters.count(MAX) != 1)
+ throw logic_error("initial uniform distribution parameter must contain max value (e.g. max= ) \n");
+ if (initialParameters.count(MIN) != 1)
+ throw logic_error("initial uniform distribution parameter must contain min value (e.g. min= ) \n");
+ float maxD = initialParameters.find(MAX)->second;
+ float minD = initialParameters.find(MIN)->second;
+ initialiseUniform(minD, maxD);
+ break;
+ }
+ case NORMAL:
+ {
+ if (initialParameters.count(MEAN) != 1)
+ throw logic_error("initial normal distribution parameter must contain mean value (e.g. mean= ) \n");
+ if (initialParameters.count(SD) != 1)
+ throw logic_error("initial normal distribution parameter must contain sdev value (e.g. sdev= ) \n");
+ float mean = initialParameters.find(MEAN)->second;
+ float sd = initialParameters.find(SD)->second;
+ initialiseNormal(mean, sd);
+ break;
+ }
+ default:
+ {
+ throw logic_error("wrong parameter value for parameter \"initialisation of dispersal trait\", must be uniform/normal \n");
+ break; //should return false
+ }
+ }
+}
+
+// ----------------------------------------------------------------------------------------
+// Dispersal initialisation options
+// ----------------------------------------------------------------------------------------
+void DispersalTrait::initialiseNormal(float mean, float sd) {
+
+ const set genePositions = pSpeciesTrait->getGenePositions();
+ short ploidy = pSpeciesTrait->getPloidy();
+
+ for (auto position : genePositions) {
+ vector> newAllelePair;
+ for (int i = 0; i < ploidy; i++) {
+ float alleleVal = pRandom->Normal(mean, sd);
+ newAllelePair.emplace_back(make_shared(alleleVal, dispDominanceFactor));
+ }
+ genes.insert(make_pair(position, newAllelePair));
+ }
+}
+
+void DispersalTrait::initialiseUniform(float min, float max) {
+
+ const set genePositions = pSpeciesTrait->getGenePositions();
+ short ploidy = pSpeciesTrait->getPloidy();
+
+ for (auto position : genePositions) {
+ vector> newAllelePair;
+ for (int i = 0; i < ploidy; i++) {
+ float alleleVal = pRandom->FRandom(min, max);
+ newAllelePair.emplace_back(make_shared(alleleVal, dispDominanceFactor));
+ }
+ genes.insert(make_pair(position, newAllelePair));
+ }
+}
+
+// ----------------------------------------------------------------------------------------
+// Dispersal gene expression options
+// ----------------------------------------------------------------------------------------
+float DispersalTrait::expressAdditive() {
+
+ float phenotype = 0.0;
+
+ for (auto const& [locus, allelePair] : genes)
+ {
+ for (const std::shared_ptr m : allelePair)
+ phenotype += m->getAlleleValue();
+ }
+ trimPhenotype(phenotype);
+ return phenotype;
+}
+
+float DispersalTrait::expressAverage() {
+
+ int positionsSize = pSpeciesTrait->getPositionsSize();
+ short ploidy = pSpeciesTrait->getPloidy();
+ float phenotype = 0.0;
+
+ for (auto const& [locus, allelePair] : genes)
+ {
+ for (auto& m : allelePair)
+ phenotype += m->getAlleleValue();
+ }
+ phenotype /= positionsSize * ploidy;
+ trimPhenotype(phenotype);
+ return phenotype;
+}
+
+void DispersalTrait::trimPhenotype(float& val) {
+ const float minPositiveVal = 1e-06;
+ switch (pSpeciesTrait->getTraitType())
+ {
+ // Values bound between 0 and 1
+ case E_D0_F: case E_D0_M: case E_D0:
+ case S_S0_F: case S_S0_M: case S_S0:
+ case KERNEL_PROBABILITY_F: case KERNEL_PROBABILITY_M: case KERNEL_PROBABILITY:
+ case CRW_STEPCORRELATION:
+ {
+ if (val < 0.0) val = 0;
+ else if (val > 1.0) val = 1.0;
+ break;
+ }
+ // Positive values
+ case KERNEL_MEANDIST_1_F: case KERNEL_MEANDIST_1_M: case KERNEL_MEANDIST_1:
+ case KERNEL_MEANDIST_2_F: case KERNEL_MEANDIST_2_M: case KERNEL_MEANDIST_2:
+ case CRW_STEPLENGTH:
+ {
+ if (val < 0.0) val = 0;
+ break;
+ }
+ // Strictly positive values
+ case E_ALPHA_F: case E_ALPHA_M: case E_ALPHA:
+ case S_ALPHA_F: case S_ALPHA_M: case S_ALPHA:
+ case SMS_ALPHADB:
+ {
+ if (val <= 0.0) val = minPositiveVal;
+ break;
+ }
+ // Minimum 1
+ case SMS_DP:
+ case SMS_GB:
+ {
+ if (val <= 1.0) val = 1.0;
+ break;
+ }
+ // Not bound
+ case E_BETA_F: case E_BETA_M: case E_BETA:
+ case S_BETA_F: case S_BETA_M: case S_BETA:
+ case SMS_BETADB:
+ {
+ break;
+ }
+ default:
+ break;
+ }
+}
+
+// ----------------------------------------------------------------------------------------
+// Get allele value at locus
+// ----------------------------------------------------------------------------------------
+float DispersalTrait::getAlleleValueAtLocus(short whichChromosome, int position) const {
+
+ auto it = genes.find(position);
+ if (it == genes.end())
+ throw runtime_error("The Dispersal locus queried for its allele value does not exist.");
+ return it->second[whichChromosome].get()->getAlleleValue();
+}
+
+float DispersalTrait::getDomCoefAtLocus(short whichChromosome, int position) const {
+ auto it = genes.find(position);
+ if (it == genes.end())
+ throw runtime_error("The genetic load locus queried for its dominance coefficient does not exist.");
+ return it->second[whichChromosome]->getDominanceCoef();
+}
+
+#ifdef UNIT_TESTS
+
+// Create a default set of alleles for testing
+//
+// Shorthand function to manually set genotypes for dispersal and
+// genetic fitness traits, instead of having to manipulate mutations.
+map>> createTestGenotype(
+ const int genomeSz, const bool isDiploid,
+ const float valAlleleA,
+ const float valAlleleB,
+ const float domCoeffA,
+ const float domCoeffB
+) {
+ vector> gene(isDiploid ? 2 : 1);
+ if (isDiploid) {
+ gene[0] = make_shared(valAlleleA, domCoeffA);
+ gene[1] = make_shared(valAlleleB, domCoeffB);
+ }
+ else {
+ gene[0] = make_shared(valAlleleA, domCoeffA);
+ }
+ map>> genotype;
+ for (int i = 0; i < genomeSz; i++) {
+ genotype.emplace(i, gene);
+ }
+ return genotype;
+}
+#endif // UNIT_TESTS
diff --git a/DispersalTrait.h b/DispersalTrait.h
new file mode 100644
index 0000000..fc5a4cc
--- /dev/null
+++ b/DispersalTrait.h
@@ -0,0 +1,129 @@
+/*----------------------------------------------------------------------------
+ *
+ * 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 .
+ *
+ * File Created by Roslyn Henry March 2023. Code adapted from NEMO (https://nemo2.sourceforge.io/)
+ --------------------------------------------------------------------------*/
+
+#ifndef DISPTRAITH
+#define DISPTRAITH
+
+#include
+#include
+#include
+#include
+
+#include "QuantitativeTrait.h"
+
+using namespace std;
+
+// Dispersal trait
+//
+// That is, all evolvable that control emigration, transfer and settlement
+class DispersalTrait : public QuantitativeTrait {
+
+public:
+
+ // Initialisation constructor, set initial values and immutable features
+ DispersalTrait(SpeciesTrait* P);
+
+ // Inheritance constructor, copies pointers to immutable features when cloning from parent
+ DispersalTrait(const DispersalTrait& T);
+
+ // Make a shallow copy to pass to offspring trait
+ // Return new pointer to new trait created by inheritance c'tor
+ // This avoids copying shared attributes: distributions and parameters
+ virtual unique_ptr clone() const override { return std::make_unique(*this); }
+
+ virtual ~DispersalTrait() { }
+
+ // Getters
+ int getNLoci() const override { return pSpeciesTrait->getPositionsSize(); }
+ float getMutationRate() const override { return pSpeciesTrait->getMutationRate(); }
+ bool isInherited() const override { return pSpeciesTrait->isInherited(); }
+
+ map>>& getGenes() { return genes; } // returning reference, receiver must be const
+
+ void mutate() override { (this->*_mutate_func_ptr) (); }
+ float express() override { return (this->*_express_func_ptr) (); }
+ void inheritGenes(const bool& fromMother, QuantitativeTrait* parent, set const& recomPositions, int startingChromosome) override;
+
+ float getAlleleValueAtLocus(short chromosome, int i) const override;
+ float getDomCoefAtLocus(short chromosome, int position) const override;
+
+#ifdef UNIT_TESTS // for testing only
+ void overwriteGenes(map>> genSeq) {
+ genes = genSeq;
+ }
+ void triggerInherit(
+ // inheritGenes requires passing a QuantitativeTrait, unfeasible in tests
+ const bool& fromMother,
+ map>> const& parentGenes,
+ set const& recomPositions,
+ int startChr) {
+ (this->*_inherit_func_ptr)(fromMother, parentGenes, recomPositions, startChr);
+ }
+#endif
+
+private:
+
+ const double dispDominanceFactor = 1.0; // no dominance for Dispersal traits (yet?)
+
+ // >
+ map>> genes;
+
+ // Initialisation
+ void initialiseUniform(float min, float max);
+ void initialiseNormal(float mean, float sd);
+
+ // Immutable features, set at initialisation
+ // and passed down to every subsequent trait copy
+ //// Species-level trait attributes, invariant across individuals
+ SpeciesTrait* pSpeciesTrait;
+ //// Species-level trait functions
+ void (DispersalTrait::* _mutate_func_ptr) (void);
+ void (DispersalTrait::* _inherit_func_ptr) (const bool& fromMother, map>> const& parent, set const& recomPositions, int parentChromosome);
+ float (DispersalTrait::* _express_func_ptr) (void);
+
+ // Possible values for immutable functions
+ //// Inheritance
+ void inheritDiploid(const bool& fromMother, map>> const& parent, set const& recomPositions, int parentChromosome);
+ void inheritHaploid(const bool& fromMother, map>> const& parent, set const& recomPositions, int parentChromosome);
+ void reInitialiseGenes(const bool& fromMother, map>> const& parentMutations, set const& recomPositions, int parentChromosome);
+ //// Mutation
+ void mutateUniform();
+ void mutateNormal();
+ void trimPhenotype(float& phenotype);
+ //// Gene expression
+ float expressAverage();
+ float expressAdditive();
+};
+
+#ifdef UNIT_TESTS
+// Test utilities
+
+map>> createTestGenotype(
+ const int genomeSz, const bool isDiploid,
+ const float valAlleleA,
+ const float valAlleleB = -99.9, // allow no value for haploids
+ const float domCoeffA = 1.0, // default for dispersal traits
+ const float domCoeffB = 1.0
+);
+#endif // UNIT_TESTS
+
+#endif // DISPTRAITH
diff --git a/FractalGenerator.cpp b/FractalGenerator.cpp
index 7a5ccc3..8650ef4 100644
--- a/FractalGenerator.cpp
+++ b/FractalGenerator.cpp
@@ -1,26 +1,26 @@
/*----------------------------------------------------------------------------
- *
- * Copyright (C) 2020 Greta Bocedi, Stephen C.F. Palmer, Justin M.J. Travis, Anne-Kathleen Malchow, Damaris Zurell
- *
+ *
+ * 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 .
- *
--------------------------------------------------------------------------*/
-
-
-//---------------------------------------------------------------------------
+
+
+
+ //---------------------------------------------------------------------------
#include "FractalGenerator.h"
//---------------------------------------------------------------------------
@@ -29,211 +29,195 @@ vector patches;
//----- Landscape creation --------------------------------------------------
-land::land(): x_coord(0), y_coord(0), value(0.0), avail(0) {}
+land::land() : x_coord(0), y_coord(0), value(0.0), avail(0) {}
bool compare(const land& z, const land& zz) //compares only the values of the cells
{
-return z.value < zz.value;
+ return z.value < zz.value;
}
-vector& fractal_landscape(int X,int Y,double Hurst,double prop,
- double maxValue,double minValue)
+vector& fractal_landscape(int X, int Y, double Hurst, double prop,
+ double maxValue, double minValue)
{
-#if RSDEBUG
-DEBUGLOG << "fractal_landscape(): X=" << X << " Y=" << Y
- << " Hurst=" << Hurst << " prop=" << prop
- << " maxValue=" << maxValue << " minValue=" << minValue
- << endl;
-#endif
-
-int ii, jj, x, y;
-int ix, iy;
-//int x0, y0, size, kx, kx2, ky, ky2;
-int kx,kx2,ky,ky2;
-
-double range; //range to draw random numbers at each iteration
-double nx, ny;
-double i, j;
-int Nx = X;
-int Ny = Y;
-
-double ran[5]; // to store each time the 5 random numbers for the random displacement
-
-int Nno; // number of cells NON suitable as habitat
-
-// exponents used to obtain the landscape dimensions
-double pow2x = log(((double)X-1.0))/log(2.0);
-double pow2y = log(((double)Y-1.0))/log(2.0);
-
-double **arena = new double *[X];
-for(ii = 0; ii < X; ii++) {
- arena[ii] = new double[Y];
-}
-patches.clear();
-// initialise all the landscape with zeroes
-for (jj = 0; jj < X; jj++) {
- for (ii = 0; ii < Y; ii++) {
- arena[jj][ii]=0;
- }
-}
+ int ii, jj, x, y;
+ int ix, iy;
+ //int x0, y0, size, kx, kx2, ky, ky2;
+ int kx, kx2, ky, ky2;
-// initialisation of the four corners
-arena[0][0] = 1.0 + pRandom->Random() * (maxValue-1.0);
-arena[0][Y-1] = 1.0 + pRandom->Random() * (maxValue-1.0);
-arena[X-1][0] = 1.0 + pRandom->Random() * (maxValue-1.0);
-arena[X-1][Y-1] = 1.0 + pRandom->Random() * (maxValue-1.0);
+ double range; //range to draw random numbers at each iteration
+ double nx, ny;
+ double i, j;
+ int Nx = X;
+ int Ny = Y;
-/////////////MIDPOINT DISPLACEMENT ALGORITHM//////////////////////////////////
-kx = (Nx-1) / 2;
-kx2 = 2 * kx;
-ky = (Ny-1) / 2;
-ky2 = 2 * ky;
+ double ran[5]; // to store each time the 5 random numbers for the random displacement
-for (ii = 0; ii < 5; ii++) //random displacement
-{
- ran[ii] = 1.0 + pRandom->Random() * (maxValue-1.0);
-}
+ int Nno; // number of cells NON suitable as habitat
-//The diamond step:
-arena[kx][ky] = ((arena[0][0] + arena[0][ky2] + arena[kx2][0] + arena[kx2][ky2])/4) + ran[0];
+ // exponents used to obtain the landscape dimensions
+ double pow2x = log(((double)X - 1.0)) / log(2.0);
+ double pow2y = log(((double)Y - 1.0)) / log(2.0);
-//The square step:
-//left
-arena[0][ky] = ((arena[0][0] +arena[0][ky2] + arena[kx][ky]) / 3) + ran[1];
-//top
-arena[kx][0] = ((arena[0][0] + arena[kx][ky] + arena[kx2][0]) / 3) + ran[2];
-//right
-arena[kx2][ky] = ((arena[kx2][0] + arena[kx][ky] + arena[kx2][ky2]) / 3) + ran[3];
-//bottom
-arena[kx][ky2] = ((arena[0][ky2] + arena[kx][ky] +arena[kx2][ky2]) / 3) + ran[4];
+ double** arena = new double* [X];
+ for (ii = 0; ii < X; ii++) {
+ arena[ii] = new double[Y];
+ }
-range = maxValue*pow(2,-Hurst);
+ patches.clear();
+ // initialise all the landscape with zeroes
+ for (jj = 0; jj < X; jj++) {
+ for (ii = 0; ii < Y; ii++) {
+ arena[jj][ii] = 0;
+ }
+ }
-i = pow2x-1;
-j = pow2y-1;
+ // initialisation of the four corners
+ arena[0][0] = 1.0 + pRandom->Random() * (maxValue - 1.0);
+ arena[0][Y - 1] = 1.0 + pRandom->Random() * (maxValue - 1.0);
+ arena[X - 1][0] = 1.0 + pRandom->Random() * (maxValue - 1.0);
+ arena[X - 1][Y - 1] = 1.0 + pRandom->Random() * (maxValue - 1.0);
-while (i > 0) {
- nx = pow(2,i)+1;
- kx = (int)((nx-1) / 2);
+ /////////////MIDPOINT DISPLACEMENT ALGORITHM//////////////////////////////////
+ kx = (Nx - 1) / 2;
kx2 = 2 * kx;
-
- ny = pow(2,j)+1;
- ky = (int)((ny-1) / 2);
+ ky = (Ny - 1) / 2;
ky2 = 2 * ky;
- ix = 0;
- while (ix <= (Nx-nx)) {
- iy = 0;
- while (iy <= (Ny-ny)) {
- for (ii = 0; ii < 5; ii++) //random displacement
- {
- ran[ii] = (int)(pRandom->Random() * 2.0 * range - range);
+ for (ii = 0; ii < 5; ii++) //random displacement
+ {
+ ran[ii] = 1.0 + pRandom->Random() * (maxValue - 1.0);
+ }
+
+ //The diamond step:
+ arena[kx][ky] = ((arena[0][0] + arena[0][ky2] + arena[kx2][0] + arena[kx2][ky2]) / 4) + ran[0];
+
+ //The square step:
+ //left
+ arena[0][ky] = ((arena[0][0] + arena[0][ky2] + arena[kx][ky]) / 3) + ran[1];
+ //top
+ arena[kx][0] = ((arena[0][0] + arena[kx][ky] + arena[kx2][0]) / 3) + ran[2];
+ //right
+ arena[kx2][ky] = ((arena[kx2][0] + arena[kx][ky] + arena[kx2][ky2]) / 3) + ran[3];
+ //bottom
+ arena[kx][ky2] = ((arena[0][ky2] + arena[kx][ky] + arena[kx2][ky2]) / 3) + ran[4];
+
+ range = maxValue * pow(2, -Hurst);
+
+ i = pow2x - 1;
+ j = pow2y - 1;
+
+ while (i > 0) {
+ nx = pow(2, i) + 1;
+ kx = (int)((nx - 1) / 2);
+ kx2 = 2 * kx;
+
+ ny = pow(2, j) + 1;
+ ky = (int)((ny - 1) / 2);
+ ky2 = 2 * ky;
+
+ ix = 0;
+ while (ix <= (Nx - nx)) {
+ iy = 0;
+ while (iy <= (Ny - ny)) {
+ for (ii = 0; ii < 5; ii++) //random displacement
+ {
+ ran[ii] = (int)(pRandom->Random() * 2.0 * range - range);
+ }
+ //The diamond step:
+
+ arena[ix + kx][iy + ky] = ((arena[ix][iy] + arena[ix][iy + ky2] + arena[ix + ky2][iy]
+ + arena[ix + kx2][iy + ky2]) / 4) + ran[0];
+ if (arena[ix + kx][iy + ky] < 1) arena[ix + kx][iy + ky] = 1;
+
+ //The square step:
+ //left
+ arena[ix][iy + ky] = ((arena[ix][iy] + arena[ix][iy + ky2] + arena[ix + kx][iy + ky]) / 3)
+ + ran[1];
+ if (arena[ix][iy + ky] < 1) arena[ix][iy + ky] = 1;
+ //top
+ arena[ix + kx][iy] = ((arena[ix][iy] + arena[ix + kx][iy + ky] + arena[ix + kx2][iy]) / 3)
+ + ran[2];
+ if (arena[ix + kx][iy] < 1) arena[ix + kx][iy] = 1;
+ //right
+ arena[ix + kx2][iy + ky] = ((arena[ix + kx2][iy] + arena[ix + kx][iy + ky] +
+ arena[ix + kx2][iy + ky2]) / 3) + ran[3];
+ if (arena[ix + kx2][iy + ky] < 1) arena[ix + kx2][iy + ky] = 1;
+ //bottom
+ arena[ix + kx][iy + ky2] = ((arena[ix][iy + ky2] + arena[ix + kx][iy + ky] +
+ arena[ix + kx2][iy + ky2]) / 3) + ran[4];
+ if (arena[ix + kx][iy + ky2] < 1) arena[ix + kx][iy + ky2] = 1;
+
+ iy += ((int)ny - 1);
}
- //The diamond step:
-
- arena[ix+kx][iy+ky] = ((arena[ix][iy] + arena[ix][iy+ky2] + arena[ix+ky2][iy]
- + arena[ix+kx2][iy+ky2])/ 4) + ran[0];
- if (arena[ix+kx][iy+ky] < 1) arena[ix+kx][iy+ky] = 1;
-
- //The square step:
- //left
- arena[ix][iy+ky] =((arena[ix][iy] +arena[ix][iy+ky2] + arena[ix+kx][iy+ky])/3)
- + ran[1];
- if (arena[ix][iy+ky] < 1) arena[ix][iy+ky] = 1;
- //top
- arena[ix+kx][iy] =((arena[ix][iy] + arena[ix+kx][iy+ky] + arena[ix+kx2][iy])/3)
- + ran[2];
- if (arena[ix+kx][iy] < 1) arena[ix+kx][iy] = 1;
- //right
- arena[ix+kx2][iy+ky] = ((arena[ix+kx2][iy] + arena[ix+kx][iy+ky] +
- arena[ix+kx2][iy+ky2]) / 3) + ran[3];
- if (arena[ix+kx2][iy+ky] < 1) arena[ix+kx2][iy+ky] = 1;
- //bottom
- arena[ix+kx][iy+ky2] = ((arena[ix][iy+ky2] + arena[ix+kx][iy+ky] +
- arena[ix+kx2][iy+ky2]) / 3) + ran[4];
- if (arena[ix+kx][iy+ky2] < 1) arena[ix+kx][iy+ky2] = 1;
-
- iy += ((int)ny-1);
+ ix += ((int)nx - 1);
}
- ix += ((int)nx-1);
- }
- if (i==j) j--;
- i--;
+ if (i == j) j--;
+ i--;
- range = range*pow(2,-Hurst); //reduce the random number range
-}
+ range = range * pow(2, -Hurst); //reduce the random number range
+ }
-// Now all the cells will be sorted and the Nno cells with the lower carrying
-// capacity will be set as matrix, i.e. with K = 0
+ // Now all the cells will be sorted and the Nno cells with the lower carrying
+ // capacity will be set as matrix, i.e. with K = 0
-land *patch;
+ land* patch;
-for (x = 0; x < X; x++) // put all the cells with their values in a vector
-{
- for (y = 0; y < Y; y++)
+ for (x = 0; x < X; x++) // put all the cells with their values in a vector
{
- patch = new land;
- patch->x_coord = x;
- patch->y_coord = y;
- patch->value = (float)arena[x][y];
- patch->avail = 1;
+ for (y = 0; y < Y; y++)
+ {
+ patch = new land;
+ patch->x_coord = x;
+ patch->y_coord = y;
+ patch->value = (float)arena[x][y];
+ patch->avail = 1;
- patches.push_back(*patch);
+ patches.push_back(*patch);
- delete patch;
+ delete patch;
+ }
}
-}
-sort(patches.begin(),patches.end(),compare); // sorts the vector
+ sort(patches.begin(), patches.end(), compare); // sorts the vector
-Nno = (int)(prop*X*Y);
-for (ii = 0; ii < Nno; ii++)
-{
- patches[ii].value = 0.0;
- patches[ii].avail = 0;
-}
+ Nno = (int)(prop * X * Y);
+ for (ii = 0; ii < Nno; ii++)
+ {
+ patches[ii].value = 0.0;
+ patches[ii].avail = 0;
+ }
-double min = (double)patches[Nno].value; // variables for the rescaling
-double max = (double)patches[X*Y-1].value;
+ double min = (double)patches[Nno].value; // variables for the rescaling
+ double max = (double)patches[X * Y - 1].value;
-double diff = max - min;
-double diffK = maxValue-minValue;
-double new_value;
+ double diff = max - min;
+ double diffK = maxValue - minValue;
+ double new_value;
-vector::iterator iter = patches.begin();
-while (iter != patches.end())
-{
- if (iter->value > 0) // rescale to a range of K between Kmin and Kmax
+ vector::iterator iter = patches.begin();
+ while (iter != patches.end())
{
- new_value = maxValue - diffK * (max - (double)iter->value) / diff;
+ if (iter->value > 0) // rescale to a range of K between Kmin and Kmax
+ {
+ new_value = maxValue - diffK * (max - (double)iter->value) / diff;
- iter->value = (float)new_value;
- }
- else iter->value = 0;
+ iter->value = (float)new_value;
+ }
+ else iter->value = 0;
- iter++;
-}
+ iter++;
+ }
-if (arena != NULL) {
-#if RSDEBUG
-//DebugGUI(("fractal_landscape(): arena=" + Int2Str((int)arena)
-// + " X=" + Int2Str(X) + " Y=" + Int2Str(Y)
-// ).c_str());
-#endif
- for(ii = 0; ii < X; ii++) {
-#if RSDEBUG
-//DebugGUI(("fractal_landscape(): ii=" + Int2Str(ii)
-// + " arena[ii]=" + Int2Str((int)arena[ii])
-// ).c_str());
-#endif
- delete[] arena[ii];
+ if (arena != NULL) {
+ for (ii = 0; ii < X; ii++) {
+ delete[] arena[ii];
+ }
+ delete[] arena;
}
- delete[] arena;
-}
-return patches;
+ return patches;
}
diff --git a/FractalGenerator.h b/FractalGenerator.h
index 24acbc7..3a2e216 100644
--- a/FractalGenerator.h
+++ b/FractalGenerator.h
@@ -1,66 +1,66 @@
/*----------------------------------------------------------------------------
- *
- * Copyright (C) 2020 Greta Bocedi, Stephen C.F. Palmer, Justin M.J. Travis, Anne-Kathleen Malchow, Damaris Zurell
- *
+ *
+ * Copyright (C) 2026 Greta Bocedi, Stephen C.F. Palmer, Justin M.J. Travis, Anne-Kathleen Malchow, Roslyn Henry, Tho 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 .
- *
--------------------------------------------------------------------------*/
-
-
-/*------------------------------------------------------------------------------
-RangeShifter v2.0 FractalGenerator
-Implements the midpoint displacement algorithm for generating a fractal Landscape,
-following:
-Saupe, D. (1988). Algorithms for random fractals. In: The Science of Fractal Images
-(eds. Pietgen, H.O. & Saupe, D.). Springer, New York, pp. 71113.
+ /*------------------------------------------------------------------------------
+ RangeShifter v2.0 FractalGenerator
-For full details of RangeShifter, please see:
-Bocedi G., Palmer S.C.F., Peer G., Heikkinen R.K., Matsinos Y.G., Watts K.
-and Travis J.M.J. (2014). RangeShifter: a platform for modelling spatial
-eco-evolutionary dynamics and species responses to environmental changes.
-Methods in Ecology and Evolution, 5, 388-396. doi: 10.1111/2041-210X.12162
+ Implements the midpoint displacement algorithm for generating a fractal Landscape,
+ following:
-Authors: Greta Bocedi & Steve Palmer, University of Aberdeen
+ Saupe, D. (1988). Algorithms for random fractals. In: The Science of Fractal Images
+ (eds. Pietgen, H.O. & Saupe, D.). Springer, New York, pp. 71113.
-Last updated: 15 July 2021 by Anne-Kathleen Malchow
-------------------------------------------------------------------------------*/
+ For full details of RangeShifter, please see:
+ Bocedi G., Palmer S.C.F., Peer G., Heikkinen R.K., Matsinos Y.G., Watts K.
+ and Travis J.M.J. (2014). RangeShifter: a platform for modelling spatial
+ eco-evolutionary dynamics and species responses to environmental changes.
+ Methods in Ecology and Evolution, 5, 388-396. doi: 10.1111/2041-210X.12162
+
+ Authors: Greta Bocedi & Steve Palmer, University of Aberdeen
+
+ Last updated: 15 July 2021 by Anne-Kathleen Malchow
+
+ ------------------------------------------------------------------------------*/
#ifndef FractalGeneratorH
#define FractalGeneratorH
#include
#include
-//using namespace std;
+ //using namespace std;
#include "Parameters.h"
class land
{
- public:
+public:
land();
int x_coord;
int y_coord;
float value;
int avail; // if 0 the patch is not available as habitat, if 1 it is
- private:
+private:
};
// IMPORTANT NOTE: X AND Y ARE TRANSPOSED, i.e. X IS THE VERTICAL CO-ORDINATE
@@ -76,10 +76,7 @@ vector& fractal_landscape(
);
bool compare(const land&, const land&);
-extern RSrandom *pRandom;
-#if RSDEBUG
-extern void DebugGUI(string);
-#endif
+extern RSrandom* pRandom;
//---------------------------------------------------------------------------
#endif
diff --git a/GeneticFitnessTrait.cpp b/GeneticFitnessTrait.cpp
new file mode 100644
index 0000000..fa1bb55
--- /dev/null
+++ b/GeneticFitnessTrait.cpp
@@ -0,0 +1,558 @@
+/*----------------------------------------------------------------------------
+ *
+ * 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 .
+ *
+ * File Created by Roslyn Henry March 2023. Code adapted from NEMO (https://nemo2.sourceforge.io/)
+ --------------------------------------------------------------------------*/
+
+#include "GeneticFitnessTrait.h"
+
+// ----------------------------------------------------------------------------------------
+// Initialisation constructor
+// Called when initialising community
+// Sets up initial values, and immutable attributes (distributions and parameters)
+// that are defined at the species-level
+// ----------------------------------------------------------------------------------------
+GeneticFitnessTrait::GeneticFitnessTrait(SpeciesTrait* P)
+{
+ pSpeciesTrait = P;
+ ExpressionType expressionType = pSpeciesTrait->getExpressionType();
+
+ _inherit_func_ptr = (pSpeciesTrait->getPloidy() == 1) ? &GeneticFitnessTrait::inheritHaploid : &GeneticFitnessTrait::inheritDiploid; //this could be changed if we wanted some alternative form of inheritance
+
+ // Set initialisation parameters
+ DistributionType initialDistribution = pSpeciesTrait->getInitialDistribution();
+ map initialParameters = pSpeciesTrait->getInitialParameters();
+ switch (initialDistribution) {
+ case UNIFORM:
+ {
+ if (initialParameters.count(MAX) != 1)
+ throw logic_error("initial uniform distribution parameter must contain max value (e.g. max= ) \n");
+ if (initialParameters.count(MIN) != 1)
+ throw logic_error("initial uniform distribution parameter must contain min value (e.g. min= ) \n");
+ break;
+ }
+ case NORMAL:
+ {
+ if (initialParameters.count(MEAN) != 1)
+ throw logic_error("initial normal distribution parameter must contain mean value (e.g. mean= ) \n");
+ if (initialParameters.count(SD) != 1)
+ throw logic_error("initial normal distribution parameter must contain sdev value (e.g. sdev= ) \n");
+ break;
+ }
+ case GAMMA:
+ {
+ if (initialParameters.count(SHAPE) != 1)
+ throw logic_error("genetic load dominance distribution set to gamma so parameters must contain one shape value (e.g. shape= ) \n");
+ if (initialParameters.count(SCALE) != 1)
+ throw logic_error("genetic load dominance distribution set to gamma so parameters must contain one scale value (e.g. scale= ) \n");
+ break;
+ }
+ case NEGEXP:
+ {
+ if (initialParameters.count(MEAN) != 1)
+ throw logic_error("genetic load dominance distribution set to negative exponential (negative decay) so parameters must contain mean value (e.g. mean= ) \n");
+ break;
+ }
+ case NONE: // initialise with default (i.e. zero) values
+ break;
+ default:
+ {
+ throw logic_error("wrong parameter value for parameter \"initialisation of dispersal traits\", must be uniform/normal \n");
+ break;
+ }
+ }
+
+ DistributionType initDomDistribution = pSpeciesTrait->getInitDomDistribution();
+ map initDomParameters = pSpeciesTrait->getInitDomParameters();
+ switch (initDomDistribution) {
+ case UNIFORM:
+ {
+ if (initDomParameters.count(MAX) != 1)
+ throw logic_error("genetic load dominance uniform distribution parameter must contain one max value (e.g. max= ) \n");
+ if (initDomParameters.count(MIN) != 1)
+ throw logic_error("genetic load dominance uniform distribution parameter must contain one min value (e.g. min= ) \n");
+ break;
+ }
+ case NORMAL:
+ {
+ if (initDomParameters.count(MEAN) != 1)
+ throw logic_error("genetic load dominance distribution set to normal so parameters must contain one mean value (e.g. mean= ) \n");
+ if (initDomParameters.count(SD) != 1)
+ throw logic_error("genetic load dominance distribution set to normal so parameters must contain one sdev value (e.g. sdev= ) \n");
+ break;
+ }
+ case GAMMA:
+ {
+ if (initDomParameters.count(SHAPE) != 1)
+ throw logic_error("genetic load dominance distribution set to gamma so parameters must contain one shape value (e.g. shape= ) \n");
+ if (initDomParameters.count(SCALE) != 1)
+ throw logic_error("genetic load dominance distribution set to gamma so parameters must contain one scale value (e.g. scale= ) \n");
+ break;
+ }
+ case NEGEXP:
+ {
+ if (initDomParameters.count(MEAN) != 1)
+ throw logic_error("genetic load dominance distribution set to negative exponential (negative decay) so parameters must contain mean value (e.g. mean= ) \n");
+ break;
+ }
+ case SCALED:
+ {
+ if (initDomParameters.count(MEAN) != 1)
+ throw logic_error("genetic load dominance distribution set to scaled, so parameters must contain mean dominance value (e.g. mean= ) \n");
+ // Set for drawing initial values
+ setScaledCoeff(initialDistribution, initialParameters);
+ break;
+ }
+ case NONE: // default values, zero-dominance coefficients
+ break;
+ default:
+ {
+ throw logic_error("wrong parameter value for genetic load dominance model, must be uniform/normal/gamma/negExp/scaled \n");
+ break;
+ }
+ }
+
+ // Draw initial values
+ initialise();
+
+ DistributionType mutationDistribution = pSpeciesTrait->getMutationDistribution();
+ map mutationParameters = pSpeciesTrait->getMutationParameters();
+ switch (mutationDistribution) {
+ case UNIFORM:
+ {
+ if (mutationParameters.count(MAX) != 1)
+ throw logic_error("genetic load mutation uniform distribution parameter must contain one max value (e.g. max= ) \n");
+ if (mutationParameters.count(MIN) != 1)
+ throw logic_error("genetic load mutation uniform distribution parameter must contain one min value (e.g. min= ) \n");
+ break;
+ }
+ case NORMAL:
+ {
+ if (mutationParameters.count(MEAN) != 1)
+ throw logic_error("genetic load mutation distribution set to normal so parameters must contain one mean value (e.g. mean= ) \n");
+ if (mutationParameters.count(SD) != 1)
+ throw logic_error("genetic load mutation distribution set to normal so parameters must contain one sdev value (e.g. sdev= ) \n");
+ break;
+ }
+ case GAMMA:
+ {
+ if (mutationParameters.count(SHAPE) != 1)
+ throw logic_error("genetic load mutation distribution set to gamma so parameters must contain one shape value (e.g. shape= ) \n");
+ if (mutationParameters.count(SCALE) != 1)
+ throw logic_error("genetic load mutation distribution set to gamma so parameters must contain one scale value (e.g. scale= ) \n");
+ break;
+ }
+ case NEGEXP:
+ {
+ if (mutationParameters.count(MEAN) != 1)
+ throw logic_error("genetic load mutation distribution set to negative exponential (negative decay) so parameters must contain one mean value (e.g. mean= ) \n");
+ break;
+ }
+ default:
+ throw logic_error("wrong parameter value for genetic load mutation model, must be uniform/normal/gamma/negExp \n");
+ }
+
+ DistributionType dominanceDistribution = pSpeciesTrait->getDominanceDistribution();
+ map dominanceParameters = pSpeciesTrait->getDominanceParameters();
+
+ switch (dominanceDistribution) {
+ case UNIFORM:
+ {
+ if (dominanceParameters.count(MAX) != 1)
+ throw logic_error("genetic load dominance uniform distribution parameter must contain one max value (e.g. max= ) \n");
+ if (dominanceParameters.count(MIN) != 1)
+ throw logic_error("genetic load dominance uniform distribution parameter must contain one min value (e.g. min= ) \n");
+ break;
+ }
+ case NORMAL:
+ {
+ if (dominanceParameters.count(MEAN) != 1)
+ throw logic_error("genetic load dominance distribution set to normal so parameters must contain one mean value (e.g. mean= ) \n");
+ if (dominanceParameters.count(SD) != 1)
+ throw logic_error("genetic load dominance distribution set to normal so parameters must contain one sdev value (e.g. sdev= ) \n");
+ break;
+ }
+ case GAMMA:
+ {
+ if (dominanceParameters.count(SHAPE) != 1)
+ throw logic_error("genetic load dominance distribution set to gamma so parameters must contain one shape value (e.g. shape= ) \n");
+ if (dominanceParameters.count(SCALE) != 1)
+ throw logic_error("genetic load dominance distribution set to gamma so parameters must contain one scale value (e.g. scale= ) \n");
+ break;
+ }
+ case NEGEXP:
+ {
+ if (dominanceParameters.count(MEAN) != 1)
+ throw logic_error("genetic load dominance distribution set to negative exponential (negative decay) so parameters must contain mean value (e.g. mean= ) \n");
+ break;
+ }
+ case SCALED:
+ {
+ if (dominanceParameters.count(MEAN) != 1)
+ throw logic_error("genetic load dominance distribution set to scaled, so parameters must contain mean dominance value (e.g. mean= ) \n");
+
+ // Set for drawing mutations (overwrite initial value)
+ setScaledCoeff(mutationDistribution, mutationParameters);
+ break;
+ }
+ default:
+ {
+ throw logic_error("wrong parameter value for genetic load dominance model, must be uniform/normal/gamma/negExp/scaled \n");
+ break;
+ }
+ }
+}
+
+// Calculate mean selection coeff s_d for calculation of k
+void GeneticFitnessTrait::setScaledCoeff(const DistributionType& selCoeffDist, const map& selCoeffParams)
+{
+ switch (selCoeffDist)
+ {
+ case UNIFORM:
+ scaledDomMeanSelCoeff = (selCoeffParams.find(MIN)->second + selCoeffParams.find(MAX)->second) / 2;
+ break;
+ case NORMAL:
+ scaledDomMeanSelCoeff = selCoeffParams.find(MEAN)->second;
+ break;
+ case GAMMA:
+ scaledDomMeanSelCoeff = selCoeffParams.find(SHAPE)->second * selCoeffParams.find(SCALE)->second;
+ break;
+ case NEGEXP:
+ scaledDomMeanSelCoeff = 1 / selCoeffParams.find(MEAN)->second;
+ break;
+ case NONE:
+ throw logic_error("Scaled dominance distribution cannot be used with default allele distribution.");
+ default: break;
+ }
+}
+
+// ----------------------------------------------------------------------------------------
+// Inheritance constructor
+// Copies immutable features from a parent trait
+// Only called via clone()
+// ----------------------------------------------------------------------------------------
+GeneticFitnessTrait::GeneticFitnessTrait(const GeneticFitnessTrait& T) :
+ pSpeciesTrait(T.pSpeciesTrait),
+ _inherit_func_ptr(T._inherit_func_ptr),
+ scaledDomMeanSelCoeff(T.scaledDomMeanSelCoeff)
+{
+ // nothing
+}
+
+void GeneticFitnessTrait::initialise() {
+ float initSelCoeff;
+ float initDomCoeff;
+ short ploidy = pSpeciesTrait->getPloidy();
+ auto initDist = pSpeciesTrait->getInitialDistribution();
+ auto initParams = pSpeciesTrait->getInitialParameters();
+ auto initDomDist = pSpeciesTrait->getInitDomDistribution();
+ auto initDomParams = pSpeciesTrait->getInitDomParameters();
+
+ const set genePositions = pSpeciesTrait->getGenePositions();
+ const set initPositions = pSpeciesTrait->getInitPositions();
+
+ for (auto position : genePositions) {
+ vector> initialGene(ploidy);
+ for (int p = 0; p < ploidy; p++) {
+ initSelCoeff = initDomCoeff = 0.0;
+ if (initPositions.contains(position)) {
+ initSelCoeff = drawSelectionCoef(initDist, initParams);
+ initDomCoeff = drawDominance(initSelCoeff, initDomDist, initDomParams);
+ }
+ initialGene[p] = make_shared(initSelCoeff, initDomCoeff);
+ }
+ genes.insert(make_pair(position, initialGene));
+ }
+}
+
+// ----------------------------------------------------------------------------------------
+// Mutate uniform
+// ----------------------------------------------------------------------------------------
+void GeneticFitnessTrait::mutate()
+{
+ const int positionsSize = pSpeciesTrait->getPositionsSize();
+ const auto& genePositions = pSpeciesTrait->getGenePositions();
+ const short ploidy = pSpeciesTrait->getPloidy();
+ const float mutationRate = pSpeciesTrait->getMutationRate();
+ float newSelectionCoef;
+ float newDominanceCoef;
+
+ auto rng = pRandom->getRNG();
+
+ for (int p = 0; p < ploidy; p++) {
+ // Determine nb of mutations
+ unsigned int NbMut = pRandom->Binomial(positionsSize, mutationRate);
+
+ if (NbMut > 0) {
+ vector mutationPositions;
+ // Draw which positions mutate
+ sample(genePositions.begin(), genePositions.end(), std::back_inserter(mutationPositions),
+ NbMut, rng);
+
+ for (int m : mutationPositions) {
+ auto it = genes.find(m);
+ if (it == genes.end())
+ throw runtime_error("Locus sampled for mutation doesn't exist.");
+ newSelectionCoef = drawSelectionCoef(
+ pSpeciesTrait->getMutationDistribution(),
+ pSpeciesTrait->getMutationParameters()
+ );
+ newDominanceCoef = drawDominance(
+ newSelectionCoef,
+ pSpeciesTrait->getDominanceDistribution(),
+ pSpeciesTrait->getDominanceParameters()
+ );
+ it->second[p] = make_shared(newSelectionCoef, newDominanceCoef);
+ }
+ }
+ }
+}
+
+// ----------------------------------------------------------------------------------------
+// get dominance value for new mutation
+// ----------------------------------------------------------------------------------------
+float GeneticFitnessTrait::drawDominance(float selCoef, const DistributionType& domDist, const map& domParams) {
+
+ float h = 0.0;
+ switch (domDist) {
+ case UNIFORM:
+ {
+ float maxD = domParams.find(MAX)->second;
+ float minD = domParams.find(MIN)->second;
+ h = pRandom->FRandom(minD, maxD);
+ break;
+ }
+ case NORMAL:
+ {
+ const float mean = domParams.find(MEAN)->second;
+ const float sd = domParams.find(SD)->second;
+ do {
+ h = static_cast(pRandom->Normal(mean, sd));
+ } while (h <= 0.0);
+ break;
+ }
+ case GAMMA:
+ {
+ const float shape = domParams.find(SHAPE)->second;
+ const float scale = domParams.find(SCALE)->second;
+ h = static_cast(pRandom->Gamma(shape, scale));
+ break;
+ }
+ case NEGEXP:
+ {
+ const float mean = domParams.find(MEAN)->second;
+ h = static_cast(pRandom->NegExp(mean));
+ break;
+ }
+ case SCALED:
+ {
+ const float h_d = domParams.find(MEAN)->second;
+ const float k = -log(2 * h_d) / scaledDomMeanSelCoeff;
+ const float max = exp(-k * selCoef);
+ h = pRandom->FRandom(0, max);
+ break;
+ }
+ case NONE:
+ {
+ // nothing, s remains 0.0
+ break;
+ }
+ default:
+ {
+ throw logic_error("wrong parameter value for genetic load dominance model, must be uniform/normal/gamma/negExp/scaled/none \n");
+ break;
+ }
+ }
+ return h;
+}
+
+// ----------------------------------------------------------------------------------------
+// Get selection coefficient for new mutation
+//
+// Selection coefficients will usually be between 0 and 1, but may,
+// if the mutation distribution enable it, take a negative value
+// down to -1 representing the effect of beneficial mutations
+// ----------------------------------------------------------------------------------------
+float GeneticFitnessTrait::drawSelectionCoef(const DistributionType& mutationDistribution, const map& mutationParameters) {
+
+ float s = 0.0; // default selection coefficient is 0
+
+ switch (mutationDistribution) {
+ case UNIFORM:
+ {
+ float maxD = mutationParameters.find(MAX)->second;
+ float minD = mutationParameters.find(MIN)->second;
+ s = pRandom->FRandom(minD, maxD); // no check here, min and max should already be constrained to valid values
+ break;
+ }
+ case NORMAL:
+ {
+ const float mean = mutationParameters.find(MEAN)->second;
+ const float sd = mutationParameters.find(SD)->second;
+ do {
+ s = static_cast(pRandom->Normal(mean, sd));
+ } while (!pSpeciesTrait->isValidTraitVal(s));
+ break;
+ }
+ case GAMMA:
+ {
+ const float shape = mutationParameters.find(SHAPE)->second;
+ const float scale = mutationParameters.find(SCALE)->second;
+ do {
+ s = static_cast(pRandom->Gamma(shape, scale));
+ } while (!pSpeciesTrait->isValidTraitVal(s));
+ break;
+ }
+ case NEGEXP:
+ {
+ const float mean = mutationParameters.find(MEAN)->second;
+ do {
+ s = static_cast(pRandom->NegExp(mean));
+ } while (!pSpeciesTrait->isValidTraitVal(s));
+ break;
+ }
+ case NONE:
+ {
+ // nothing, s remains 0.0
+ break;
+ }
+ default:
+ {
+ throw logic_error("wrong parameter value for genetic load mutation model, must be uniform/normal/gamma/negExp/scaled/none \n");
+ break;
+ }
+ }
+ return s;
+}
+
+
+// ----------------------------------------------------------------------------------------
+// Wrapper to inheritance function
+// ----------------------------------------------------------------------------------------
+void GeneticFitnessTrait::inheritGenes(const bool& fromMother, QuantitativeTrait* parentTrait, set const& recomPositions, int startingChromosome)
+{
+ auto parentCast = dynamic_cast (parentTrait); // must convert QuantitativeTrait to GeneticFitnessTrait
+ const auto& parent_seq = parentCast->getGenes();
+ (this->*_inherit_func_ptr) (fromMother, parent_seq, recomPositions, startingChromosome);
+}
+
+// ----------------------------------------------------------------------------------------
+// Inheritance for diploid, sexual species
+// Called once for each parent. Given a list of recombinant sites,
+// populates offspring genes with appropriate parent alleles
+// Assumes mother genes are inherited first
+// ----------------------------------------------------------------------------------------
+void GeneticFitnessTrait::inheritDiploid(const bool& fromMother, map>> const& parentGenes, set const& recomPositions, int parentChromosome) {
+
+ const int lastPosition = parentGenes.rbegin()->first;
+ auto recomIt = recomPositions.lower_bound(parentGenes.begin()->first);
+ // If no recombination sites, only breakpoint is last position
+ // i.e., no recombination occurs
+ int nextBreakpoint = recomIt == recomPositions.end() ? lastPosition : *recomIt;
+
+ // Is the first parent gene position already recombinant?
+ auto distance = std::distance(recomPositions.begin(), recomIt);
+ if (distance % 2 != 0)
+ parentChromosome = 1 - parentChromosome; // switch chromosome
+
+ for (auto const& [locus, allelePair] : parentGenes) {
+
+ // Switch chromosome if locus is past recombination site
+ while (locus > nextBreakpoint) {
+ parentChromosome = 1 - parentChromosome;
+ std::advance(recomIt, 1); // go to next recombination site
+ nextBreakpoint = recomIt == recomPositions.end() ? lastPosition : *recomIt;
+ }
+
+ if (locus <= nextBreakpoint) {
+ auto& parentAllele = allelePair[parentChromosome];
+ auto itGene = genes.find(locus);
+ if (itGene == genes.end()) {
+ // locus does not exist yet, create and initialise it
+ if (!fromMother) throw runtime_error("Father-inherited locus does not exist.");
+ vector> newAllelePair(2); // always diploid
+ newAllelePair[sex_t::FEM] = parentAllele;
+ genes.insert(make_pair(locus, newAllelePair));
+ }
+ else { // father, locus already exists
+ if (fromMother) throw runtime_error("Mother-inherited locus already exists.");
+ itGene->second[sex_t::MAL] = parentAllele;
+ }
+ }
+ }
+}
+
+// ----------------------------------------------------------------------------------------
+// Inheritance for haploid, asexual species
+// Simply pass down parent genes
+// Arguments are still needed to match overloaded function in base class
+// ----------------------------------------------------------------------------------------
+void GeneticFitnessTrait::inheritHaploid(const bool& fromMother, map>> const& parentGenes, set const& recomPositions, int parentChromosome)
+{
+ genes = parentGenes;
+}
+
+// ----------------------------------------------------------------------------------------
+// Expression genetic load
+// ----------------------------------------------------------------------------------------
+float GeneticFitnessTrait::express() {
+
+ float phenotype = 1.0; // base chance of viability
+ float sA, sB, hA, hB, sumDomCoeffs, hLocus;
+
+ for (auto const& [locus, pAllelePair] : genes)
+ {
+ shared_ptr pAlleleA = pAllelePair[0] == 0 ? wildType : pAllelePair[0];
+
+ sA = pAlleleA->getAlleleValue();
+ hA = pAlleleA->getDominanceCoef();
+
+ if (pSpeciesTrait->getPloidy() == 2) {
+ shared_ptr pAlleleB = pAllelePair[1] == 0 ? wildType : pAllelePair[1];
+ sB = pAlleleB->getAlleleValue();
+ hB = pAlleleB->getDominanceCoef();
+
+ sumDomCoeffs = hA + hB;
+ hLocus = sumDomCoeffs == 0.0 ? 0.5 : hA / sumDomCoeffs;
+ phenotype *= 1 - hLocus * sA - (1 - hLocus) * sB;
+ }
+ else {
+ phenotype *= 1 - sA;
+ }
+ }
+ return phenotype;
+}
+
+// ----------------------------------------------------------------------------------------
+// Get allele value at locus
+// ----------------------------------------------------------------------------------------
+float GeneticFitnessTrait::getAlleleValueAtLocus(short whichChromosome, int position) const {
+
+ auto it = genes.find(position);
+ if (it == genes.end())
+ throw runtime_error("The genetic load locus queried for its allele value does not exist.");
+ return it->second[whichChromosome] == 0 ? wildType->getAlleleValue() : it->second[whichChromosome]->getAlleleValue();
+}
+
+float GeneticFitnessTrait::getDomCoefAtLocus(short whichChromosome, int position) const {
+ auto it = genes.find(position);
+ if (it == genes.end())
+ throw runtime_error("The genetic load locus queried for its dominance coefficient does not exist.");
+ return it->second[whichChromosome] == 0 ? wildType->getDominanceCoef() : it->second[whichChromosome]->getDominanceCoef();
+}
diff --git a/GeneticFitnessTrait.h b/GeneticFitnessTrait.h
new file mode 100644
index 0000000..cb81eff
--- /dev/null
+++ b/GeneticFitnessTrait.h
@@ -0,0 +1,111 @@
+/*----------------------------------------------------------------------------
+ *
+ * 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 .
+ *
+ * File Created by Roslyn Henry March 2023. Code adapted from NEMO (https://nemo2.sourceforge.io/)
+ --------------------------------------------------------------------------*/
+
+#ifndef GENETICFITNESSH
+#define GENETICFITNESSH
+
+#include
+#include
+#include
+#include