diff --git a/src/RScore/.github/workflows/check.yml b/src/RScore/.github/workflows/check.yml new file mode 100644 index 0000000..2e86587 --- /dev/null +++ b/src/RScore/.github/workflows/check.yml @@ -0,0 +1,26 @@ +name: check +on: push + +jobs: + check: + strategy: + matrix: + os: [windows-latest, ubuntu-latest, macos-latest] + include: + - os: windows-latest + path_to_exe: ./build/Debug/RScore.exe + - os: ubuntu-latest + path_to_exe: ./build/RScore + - os: macos-latest + path_to_exe: ./build/RScore + runs-on: ${{ matrix.os }} + steps: + - uses: actions/checkout@v3 + + - name: build + run: | + mkdir build && cd build + cmake ../ && cmake --build . + + - name: run + run: ${{ matrix.path_to_exe }} diff --git a/src/RScore/Allele.h b/src/RScore/Allele.h index 7f8a0a8..1b00706 100644 --- a/src/RScore/Allele.h +++ b/src/RScore/Allele.h @@ -11,6 +11,6 @@ class Allele { ~Allele() {} float getAlleleValue() const { return value; }; float getDominanceCoef() const { return dominance; }; - float getId() const { return id; } + int getId() const { return id; } }; #endif \ No newline at end of file diff --git a/src/RScore/CMakeLists.txt b/src/RScore/CMakeLists.txt index a100830..d8fe69f 100644 --- a/src/RScore/CMakeLists.txt +++ b/src/RScore/CMakeLists.txt @@ -1,11 +1,30 @@ # Config file for compilation with CMake -# Only relevant for Batch mode, but this file needs to live in the RScore folder -add_library(RScore Species.cpp Cell.cpp Community.cpp FractalGenerator.cpp GeneticLoad.cpp Individual.cpp Landscape.cpp Model.cpp NeutralStatsManager.cpp Parameters.cpp Patch.cpp Population.cpp QTLTrait.cpp RSrandom.cpp SNPTrait.cpp SubCommunity.cpp) +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) + # specify the C++ standard + set(CMAKE_CXX_STANDARD 20) + set(CMAKE_CXX_STANDARD_REQUIRED True) + add_executable(RScore Main.cpp Species.cpp Cell.cpp Community.cpp FractalGenerator.cpp GeneticLoad.cpp Individual.cpp Landscape.cpp Model.cpp NeutralStatsManager.cpp Parameters.cpp Patch.cpp Population.cpp QTLTrait.cpp RSrandom.cpp SNPTrait.cpp SpeciesTrait.cpp SubCommunity.cpp Utils.cpp) +else() # that is, RScore compiled as library within RangeShifter_batch + add_library(RScore Species.cpp Cell.cpp Community.cpp FractalGenerator.cpp GeneticLoad.cpp Individual.cpp Landscape.cpp Model.cpp NeutralStatsManager.cpp Parameters.cpp Patch.cpp Population.cpp QTLTrait.cpp RSrandom.cpp SNPTrait.cpp SpeciesTrait.cpp SubCommunity.cpp Utils.cpp) +endif() # pass config definitions to compiler +target_compile_definitions(RScore PRIVATE RSWIN64) + +# enable LINUX_CLUSTER macro on Linux + macOS if(CMAKE_SYSTEM_NAME STREQUAL "Linux" OR CMAKE_SYSTEM_NAME STREQUAL "Darwin") - target_compile_definitions(RScore PRIVATE "RSDEBUG" "RSWIN64" "LINUX_CLUSTER") -else() # Windows - target_compile_definitions(RScore PRIVATE "RSDEBUG" "RSWIN64") + 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 diff --git a/src/RScore/CONTRIBUTING.md b/src/RScore/CONTRIBUTING.md new file mode 100644 index 0000000..4b180d3 --- /dev/null +++ b/src/RScore/CONTRIBUTING.md @@ -0,0 +1,81 @@ +# The RangeShifter platform - An eco-evolutionary modelling framework + +## How to contribute + +Thank you for your interest in contributing to the RangeShifter platform. +In this document we will give you guidance on how to contribute to the RangeShifter project regarding issues, bug fixing and adding new features. + +## Repo structure + +![Rangeshifter repo structure](RS_repos.png) + +RangeShifter is distributed with three user interfaces, each living in their own repo: + +- the RangeShifter GUI (clickable Windows interface)* +- RangeShifter Batch Mode (command line interface) +- the RangeShiftR package (R interface) + +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). + +*The RangeShifter GUI is currently being rewritten, and is not open source yet. + +## Roles + +#### 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. + +#### Developers + +Regular contributors and members of the [RangeShifter development team](https://github.com/orgs/RangeShifter/people), including maintainers. + +#### Contributors + +Anyone who whishes to make changes to RangeShifter's code, including regular developers. + +## Branching policy + +![](branches.png) + +*Check out the [Git cheatsheet](https://github.com/RangeShifter/RScore/git_cheatsheet.md) for a reminder on the main git commands* + +This policy applies to RScore and all three RangeShifter interfaces. +RangeShifter uses the following branching structure: + +- `main` is the default branch, where stable releases live. Because it contains the version of RangeShifter that users normally interact with, it must be stable and build at all times. +Only maintainers should make significant changes to `main`, normally by merging `develop` into `main` to make newly developed features available to users, and marking a release while doing so. +- `develop` is the development branch containing new, in-development features. It is the reference branch for all developers. Contributors may make small changes directly to `develop` but should ensure that new changes do not break the build. If one happens to break `develop`, it should be their top priority to fix it as this will disrupt the work of all other contributors. +Larger changes should instead be developed on feature branches. +- Larger changes should be first developed on feature (e.g. `cmake`, `mutualism`, etc.) or contributor (e.g., `theo`) branches. Contributors are welcome to experiment and break such branches at any time, as this will not impact users or other contributors. + +When progress is deemed satisfactory, changes can be brought to `develop`. Please open a pull request on GitHub, and assign at least one maintainer as a reviewer. As a pre-requisite, RangeShifter must build on the branch before merging. Please enter a descriptive title and use the description field to describe what you have changed. + +In the meantime, we encourage contributors to work in small and frequent commits, and to merge `develop` into their branch often to update their branch with newest changes. + +### Contributing to RangeShifter core code + +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 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. + +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 + +Do you have an idea of a new feature in the RangeShifter platform that should be integrated and is of use for other RangeShifter users? +Please get in touch with the RangeShifter development team (rangeshiftr@uni-potsdam.de) to discuss a collaboration. + +⚠️ We advise to contact the developer team as early as possible if you plan on implementing a new feature. This could prevent simultaneous development of the same feature and coordinate potential joint development. + +Alternatively, proceed as with the bug fix above: create your own branch or fork _from `develop`_ and work from there, and submit a pull request when your new features are ready to join the core code. +We recommend that you update your branch regularly to new changes on `develop` (using `git merge develop`) to reduce the risk of merge conflicts or your version getting out-of-touch in the late stages of development. +We also recommend that you work in small commits, as this makes the code easier to debug, and makes it easier for maintainers to understand your contributions when reviewing a pull request. diff --git a/src/RScore/Cell.cpp b/src/RScore/Cell.cpp index 2f8c1b3..fac40c8 100644 --- a/src/RScore/Cell.cpp +++ b/src/RScore/Cell.cpp @@ -59,18 +59,12 @@ Cell::Cell(int xx, int yy, intptr patch, float hab) } 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; } -#if RSDEBUG - //DEBUGLOG << "Cell::~Cell(): deleted" << endl; -#endif } void Cell::setHabIndex(short hx) { @@ -118,12 +112,6 @@ void Cell::setPatch(intptr p) { } intptr 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; } diff --git a/src/RScore/Community.cpp b/src/RScore/Community.cpp index c05e023..a8f7673 100644 --- a/src/RScore/Community.cpp +++ b/src/RScore/Community.cpp @@ -64,7 +64,7 @@ void Community::initialise(Species* pSpecies, int year) int nsubcomms, npatches, ndistcells, spratio, patchnum, rr = 0; locn distloc; patchData pch; - patchLimits limits; + patchLimits limits = patchLimits(); intptr ppatch, subcomm; std::vector subcomms; std::vector selected; @@ -257,18 +257,12 @@ 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); while (iind.year == year) { -#if RSDEBUG - //DEBUGLOG << "Community::initialise(): year=" << year - // << " indIx=" << indIx << " iind.year=" << iind.year - // << " iind.patchID=" << iind.patchID << " iind.x=" << iind.x << " iind.y=" << iind.y - // << " iind.sex=" << iind.sex << " iind.age=" << iind.age << " iind.stage=" << iind.stage - // << endl; -#endif if (ppLand.patchModel) { if (pLandscape->existsPatch(iind.patchID)) { pPatch = pLandscape->findPatch(iind.patchID); @@ -512,7 +506,7 @@ void Community::dispersal(short landIx) #endif // SEASONAL || RS_RCPP { #if RSDEBUG - int t0, t1, t2; + time_t t0, t1, t2; t0 = time(0); #endif @@ -537,21 +531,12 @@ void Community::dispersal(short landIx) for (int i = 0; i < nsubcomms; i++) { // all populations subComms[i]->resetPossSettlers(); } -#if RSDEBUG - //DEBUGLOG << "Community::dispersal() 1111: ndispersers=" << ndispersers << endl; -#endif #if RS_RCPP // included also SEASONAL ndispersers = matrix->transfer(pLandscape, landIx, nextseason); #else ndispersers = matrix->transfer(pLandscape, landIx); #endif // SEASONAL || RS_RCPP -#if RSDEBUG - //DEBUGLOG << "Community::dispersal() 2222: ndispersers=" << ndispersers << endl; -#endif matrix->completeDispersal(pLandscape, sim.outConnect); -#if RSDEBUG - //DEBUGLOG << "Community::dispersal() 3333: ndispersers=" << ndispersers << endl; -#endif } while (ndispersers > 0); #if RSDEBUG @@ -560,11 +545,6 @@ void Community::dispersal(short landIx) DEBUGLOG << "Community::dispersal(): transfer time=" << t2 - t1 << endl; #endif -#if RSDEBUG - //int t3 = time(0); - //DEBUGLOG << "Community::dispersal(): completion time = " << t3-t2 << endl; -#endif - } void Community::survival(short part, short option0, short option1) @@ -652,7 +632,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; @@ -663,25 +643,9 @@ commStats Community::getStats(void) patchPop = subComms[i]->getPopStats(); s.ninds += patchPop.nInds; s.nnonjuvs += patchPop.nNonJuvs; -#if RSDEBUG - //DEBUGLOG << "Community::getStats(): i = " << i - // << " pSpecies = " << patchPop.pSpecies << " pPatch = " << patchPop.pPatch - // << " nInds = " << patchPop.nInds << endl; -#endif if (patchPop.pPatch != 0) { // not the matrix patch -#if RSDEBUG - //DEBUGLOG << "Community::getStats(): i = " << i - // << " patchNum = " << patchPop.pPatch->getPatchNum() << endl; -#endif if (patchPop.pPatch->getPatchNum() != 0) { // not matrix patch localK = patchPop.pPatch->getK(); -#if RSDEBUG - //DEBUGLOG << "Community::getStats(): i= " << i - // << " pSpecies= " << patchPop.pSpecies << " pPatch= " << patchPop.pPatch - // << " patchNum= " << patchPop.pPatch->getPatchNum() << " localK= " << localK - // << " nInds= " << patchPop.nInds << " breeding= " << (int)patchPop.breeding - // << endl; -#endif if (localK > 0.0) s.suitable++; if (patchPop.nInds > 0 && patchPop.breeding) { s.occupied++; @@ -920,9 +884,9 @@ 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; + traitCanvas tcanv = traitCanvas(); int ngenes, popsize; tcanv.pcanvas[0] = NULL; @@ -962,17 +926,6 @@ void Community::outRange(Species* pSpecies, int rep, int yr, int gen) ts.sumS0[j] += scts.sumS0[j]; ts.ssqS0[j] += scts.ssqS0[j]; ts.sumAlphaS[j] += scts.sumAlphaS[j]; ts.ssqAlphaS[j] += scts.ssqAlphaS[j]; ts.sumBetaS[j] += scts.sumBetaS[j]; ts.ssqBetaS[j] += scts.ssqBetaS[j]; -#if RSDEBUG - //DEBUGLOG << "Community::outRange(): i=" << i << " j=" << j - // << " scts.ninds[j]=" << scts.ninds[j] - // << " scts.sumD0[j]=" << scts.sumD0[j] - // << " scts.ssqD0[j]=" << scts.ssqD0[j] - // << endl; - //DEBUGLOG << "Community::outRange(): i=" << i << " j=" << j - // << " ts.ninds[j]=" << ts.ninds[j] - // << " ts.sumD0[j]=" << ts.sumD0[j] << " ts.ssqD0[j]=" << ts.ssqD0[j] - // << endl; -#endif } } @@ -1011,16 +964,6 @@ void Community::outRange(Species* pSpecies, int rep, int yr, int gen) sdD0[g] = sdAlpha[g] = sdBeta[g] = 0.0; } } -#if RSDEBUG - //DEBUGLOG << "Community::outRange(): ngenes=" << ngenes << " g=" << g - // << " ts.ninds[g]=" << ts.ninds[g] - // << " ts.sumD0[g]=" << ts.sumD0[g] - // << " ts.ssqD0[g]=" << ts.ssqD0[g] - // << endl; - //DEBUGLOG << "Community::outRange(): popsize=" << popsize - // << " mnD0[g]" << mnD0[g] << " sdD0[g]" << sdD0[g] - // << endl; -#endif } if (emig.sexDep) { outrange << "\t" << mnD0[0] << "\t" << sdD0[0]; @@ -1098,17 +1041,6 @@ void Community::outRange(Species* pSpecies, int rep, int yr, int gen) if (sdBetaDB[g] > 0.0) sdBetaDB[g] = sqrt(sdBetaDB[g]); else sdBetaDB[g] = 0.0; } } -#if RSDEBUG - //DEBUGLOG << "Community::outRange(): ngenes=" << ngenes << " g=" << g - // << " ts.ninds[g]=" << ts.ninds[g] - // << " ts.sumDP[g]=" << ts.sumDP[g] << " ts.ssqDP[g]=" << ts.ssqDP[g] - // << " ts.sumGB[g]=" << ts.sumGB[g] << " ts.ssqGB[g]=" << ts.ssqGB[g] - // << endl; - //DEBUGLOG << "Community::outRange(): popsize=" << popsize - // << " mnDP[g]" << mnDP[g] << " sdDP[g]" << sdDP[g] - // << " mnGB[g]" << mnGB[g] << " sdGB[g]" << sdGB[g] - // << endl; -#endif } if (trfr.moveModel) { if (trfr.moveType == 1) { @@ -1282,35 +1214,17 @@ void Community::outOccSuit(bool view) { simParams sim = paramsSim->getSim(); //streamsize prec = outsuit.precision(); -#if RSDEBUG -//DEBUGLOG << "Community::outOccSuit(): sim.reps=" << sim.reps -// << " sim.years=" << sim.years << " sim.outInt=" << sim.outInt << endl; -#endif for (int i = 0; i < (sim.years / sim.outIntOcc) + 1; i++) { sum = ss = 0.0; for (int rep = 0; rep < sim.reps; rep++) { sum += occSuit[i][rep]; ss += occSuit[i][rep] * occSuit[i][rep]; -#if RSDEBUG - //DEBUGLOG << "Community::outOccSuit(): i=" << i << " rep=" << rep - // << " occSuit[i][rep]=" << occSuit[i][rep] - // << " sum=" << sum << " ss=" << ss - // << endl; -#endif } mean = sum / (double)sim.reps; sd = (ss - (sum * sum / (double)sim.reps)) / (double)(sim.reps - 1); -#if RSDEBUG - //DEBUGLOG << "Community::outOccSuit(): i=" << i - // << " mean=" << mean << " sd=" << sd << endl; -#endif if (sd > 0.0) sd = sqrt(sd); else sd = 0.0; se = sd / sqrt((double)(sim.reps)); -#if RSDEBUG - //DEBUGLOG << "Community::outOccSuit(): i=" << i - // << " sd=" << sd << " se=" << se << endl; -#endif // outsuit << i*sim.outInt << "\t" << mean << "\t" << se << endl; // if (view) viewOccSuit(i*sim.outInt,mean,se); @@ -1735,7 +1649,6 @@ Rcpp::IntegerMatrix Community::addYearToPopList(int rep, int yr) { // TODO: def intptr subcomm = 0; SubCommunity* pSubComm = 0; popStats pop; - //pop.breeding = false; pop.nInds = pop.nAdults = pop.nNonJuvs = 0; for (int y = 0; y < ppLand.dimY; y++) { @@ -1977,8 +1890,8 @@ void Community::writeWCPerLocusFstatFile(Species* pSpecies, const int yr, const const auto pPop = (Population*)patch->getPopn((intptr)pSpecies); int popSize = pPop->sampleSize(); int het = 0; - for (unsigned int a = 0; a < nAlleles; ++a) { - het += pPop->getHetero(thisLocus, a); + for (int a = 0; a < nAlleles; ++a) { + het += static_cast(pPop->getHetero(thisLocus, a)); // not sure why this returns a double } outperlocusfstat << "\t" << het / (2.0 * popSize); } diff --git a/src/RScore/Community.h b/src/RScore/Community.h index bbfde1a..0c5d3e4 100644 --- a/src/RScore/Community.h +++ b/src/RScore/Community.h @@ -107,7 +107,7 @@ class Community { // 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 + // 1 - development and survival ); void ageIncrement(void); int totalInds(void); diff --git a/src/RScore/FractalGenerator.cpp b/src/RScore/FractalGenerator.cpp index 9dda608..c0a0392 100644 --- a/src/RScore/FractalGenerator.cpp +++ b/src/RScore/FractalGenerator.cpp @@ -211,17 +211,7 @@ vector& fractal_landscape(int X, int Y, double Hurst, double prop, } 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]; } delete[] arena; diff --git a/src/RScore/GeneticLoad.cpp b/src/RScore/GeneticLoad.cpp index 7e45c19..aa93e8d 100644 --- a/src/RScore/GeneticLoad.cpp +++ b/src/RScore/GeneticLoad.cpp @@ -21,39 +21,39 @@ GeneticLoad::GeneticLoad(SpeciesTrait* P) switch (mutationDistribution) { case UNIFORM: { - if (!mutationParameters.count(MAX)) - cout << endl << ("Error:: adaptive mutation uniform distribution parameter must contain max value (e.g. max= ) \n"); + if (mutationParameters.count(MAX) != 1) + cout << endl << ("Error:: adaptive mutation uniform distribution parameter must contain one max value (e.g. max= ) \n"); - if (!mutationParameters.count(MIN)) - cout << endl << ("Error:: adaptive mutation uniform distribution parameter must contain min value (e.g. min= ) \n"); + if (mutationParameters.count(MIN) != 1) + cout << endl << ("Error:: adaptive mutation uniform distribution parameter must contain one min value (e.g. min= ) \n"); break; } case NORMAL: { - if (!mutationParameters.count(MEAN)) - cout << endl << ("Error:: adaptive mutation distribution set to normal so parameters must contain mean value (e.g. mean= ) \n"); + if (mutationParameters.count(MEAN) != 1) + cout << endl << ("Error:: adaptive mutation distribution set to normal so parameters must contain one mean value (e.g. mean= ) \n"); - if (!mutationParameters.count(SDEV)) - cout << endl << ("Error:: adaptive mutation distribution set to normal so parameters must contain sdev value (e.g. sdev= ) \n"); + if (mutationParameters.count(SDEV) != 1) + cout << endl << ("Error:: adaptive mutation distribution set to normal so parameters must contain one sdev value (e.g. sdev= ) \n"); break; } case GAMMA: { - if (!mutationParameters.count(SHAPE)) - cout << endl << ("Error:: adaptive mutation distribution set to gamma so parameters must contain shape value (e.g. shape= ) \n"); + if (mutationParameters.count(SHAPE) != 1) + cout << endl << ("Error:: adaptive mutation distribution set to gamma so parameters must contain one shape value (e.g. shape= ) \n"); - if (!mutationParameters.count(SCALE)) - cout << endl << ("Error:: adaptive mutation distribution set to gamma so parameters must contain scale value (e.g. scale= ) \n"); + if (mutationParameters.count(SCALE) != 1) + cout << endl << ("Error:: adaptive mutation distribution set to gamma so parameters must contain one scale value (e.g. scale= ) \n"); break; } case NEGEXP: { - if (!mutationParameters.count(MEAN)) - cout << endl << ("Error:: adaptive mutation distribution set to negative exponential (negative decay) so parameters must contain mean value (e.g. mean= ) \n"); + if (mutationParameters.count(MEAN) != 1) + cout << endl << ("Error:: adaptive mutation distribution set to negative exponential (negative decay) so parameters must contain one mean value (e.g. mean= ) \n"); break; } @@ -70,38 +70,38 @@ GeneticLoad::GeneticLoad(SpeciesTrait* P) switch (dominanceDistribution) { case UNIFORM: { - if (!dominanceParameters.count(MAX)) - cout << endl << ("Error:: adaptive dominance uniform distribution parameter must contain max value (e.g. max= ) \n"); + if (dominanceParameters.count(MAX) != 1) + cout << endl << ("Error:: adaptive dominance uniform distribution parameter must contain one max value (e.g. max= ) \n"); - if (!dominanceParameters.count(MIN)) - cout << endl << ("Error:: adaptive dominance uniform distribution parameter must contain min value (e.g. min= ) \n"); + if (dominanceParameters.count(MIN) != 1) + cout << endl << ("Error:: adaptive dominance uniform distribution parameter must contain one min value (e.g. min= ) \n"); break; } case NORMAL: { - if (!dominanceParameters.count(MEAN)) - cout << endl << ("Error:: adaptive dominance distribution set to normal so parameters must contain mean value (e.g. mean= ) \n"); + if (dominanceParameters.count(MEAN) != 1) + cout << endl << ("Error:: adaptive dominance distribution set to normal so parameters must contain one mean value (e.g. mean= ) \n"); - if (!dominanceParameters.count(SDEV)) - cout << endl << ("Error:: adaptive dominance distribution set to normal so parameters must contain sdev value (e.g. sdev= ) \n"); + if (dominanceParameters.count(SDEV) != 1) + cout << endl << ("Error:: adaptive dominance distribution set to normal so parameters must contain one sdev value (e.g. sdev= ) \n"); break; } case GAMMA: { - if (!dominanceParameters.count(SHAPE)) - cout << endl << ("Error:: adaptive dominance distribution set to gamma so parameters must contain shape value (e.g. shape= ) \n"); + if (dominanceParameters.count(SHAPE) != 1) + cout << endl << ("Error:: adaptive dominance distribution set to gamma so parameters must contain one shape value (e.g. shape= ) \n"); - if (!dominanceParameters.count(SCALE)) - cout << endl << ("Error:: adaptive dominance distribution set to gamma so parameters must contain scale value (e.g. scale= ) \n"); + if (dominanceParameters.count(SCALE) != 1) + cout << endl << ("Error:: adaptive dominance distribution set to gamma so parameters must contain one scale value (e.g. scale= ) \n"); break; } case NEGEXP: { - if (!dominanceParameters.count(MEAN)) + if (dominanceParameters.count(MEAN) != 1) cout << endl << ("Error:: adaptive dominance distribution set to negative exponential (negative decay) so parameters must contain mean value (e.g. mean= ) \n"); break; @@ -211,27 +211,27 @@ float GeneticLoad::drawDominance(float selCoef) { { const float mean = dominanceParameters.find(MEAN)->second; const float sd = dominanceParameters.find(SDEV)->second; - h = pRandom->Normal(mean, sd); + h = static_cast(pRandom->Normal(mean, sd)); break; } case GAMMA: { const float shape = dominanceParameters.find(SHAPE)->second; const float scale = dominanceParameters.find(SCALE)->second; - h = pRandom->Gamma(shape, scale); + h = static_cast(pRandom->Gamma(shape, scale)); break; } case NEGEXP: { const float mean = dominanceParameters.find(MEAN)->second; - h = pRandom->NegExp(mean); + h = static_cast(pRandom->NegExp(mean)); break; } case SCALED: { const float min = 0; - const float max = exp((-log(2 * 0.36) / 0.05) * selCoef); - h = pRandom->FRandom(min, max); + const float max = static_cast(exp((-log(2 * 0.36) / 0.05) * selCoef)); + h = static_cast(pRandom->FRandom(min, max)); break; } @@ -270,7 +270,7 @@ float GeneticLoad::drawSelectionCoef() { { const float mean = mutationParameters.find(MEAN)->second; const float sd = mutationParameters.find(SDEV)->second; - s = pRandom->Normal(mean, sd); + s = static_cast(pRandom->Normal(mean, sd)); break; } @@ -278,13 +278,13 @@ float GeneticLoad::drawSelectionCoef() { { const float shape = mutationParameters.find(SHAPE)->second; const float scale = mutationParameters.find(SCALE)->second; - s = pRandom->Gamma(shape, scale); + s = static_cast(pRandom->Gamma(shape, scale)); break; } case NEGEXP: { const float mean = mutationParameters.find(MEAN)->second; - s = pRandom->NegExp(mean); + s = static_cast(pRandom->NegExp(mean)); break; } default: @@ -315,11 +315,11 @@ void GeneticLoad::inheritDiploid(sex_t whichChromosome, mapfirst); - unsigned int nextBreakpoint = *it; + int nextBreakpoint = *it; auto distance = std::distance(recomPositions.begin(), it); if (distance % 2 != 0) - parentChromosome = !parentChromosome; // switch to the other one + parentChromosome = 1 - parentChromosome; // switch to the other one // use 1-parentChromosome, or switch to a sex_t ? for (auto const& [locus, allelePair] : parentGenes) { @@ -327,7 +327,7 @@ void GeneticLoad::inheritDiploid(sex_t whichChromosome, map nextBreakpoint) { std::advance(it, 1); nextBreakpoint = *it; - parentChromosome = !parentChromosome; // switch to the other one + parentChromosome = 1 - parentChromosome; // switch to the other one } if (locus <= nextBreakpoint) { @@ -363,8 +363,8 @@ float GeneticLoad::express() { for (auto const& [locus, pAllelePair] : genes) { - auto pAlleleLeft = (!pAllelePair[0]) ? wildType : pAllelePair[0]; - auto pAlleleRight = (!pAllelePair[1]) ? wildType : pAllelePair[1]; + shared_ptr pAlleleLeft = (!pAllelePair[0]) ? wildType : pAllelePair[0]; + shared_ptr pAlleleRight = (!pAllelePair[1]) ? wildType : pAllelePair[1]; if (pAlleleLeft.get()->getId() != pAlleleRight.get()->getId()) // heterozygote { @@ -391,9 +391,9 @@ bool GeneticLoad::isHeterozygoteAtLocus(int locus) const { if (it == genes.end()) //not found so must be wildtype homozygous return false; else { - auto a = (!it->second[0]) ? wildType : it->second[0]; - auto b = (!it->second[1]) ? wildType : it->second[1]; - return a != b; + shared_ptr alleleRight = (!it->second[0]) ? wildType : it->second[0]; + shared_ptr alleleLeft = (!it->second[1]) ? wildType : it->second[1]; + return alleleRight != alleleLeft; } } @@ -407,8 +407,8 @@ int GeneticLoad::countHeterozygoteLoci() const { int count = 0; for (auto const& [locus, allelePair] : genes) { - auto alleleLeft = (!allelePair[0]) ? wildType : allelePair[0]; - auto alleleRight = (!allelePair[1]) ? wildType : allelePair[1]; + shared_ptr alleleLeft = (!allelePair[0]) ? wildType : allelePair[0]; + shared_ptr alleleRight = (!allelePair[1]) ? wildType : allelePair[1]; count += alleleLeft != alleleRight; } return count; diff --git a/src/RScore/Individual.cpp b/src/RScore/Individual.cpp index 25dc99f..8a9fbb8 100644 --- a/src/RScore/Individual.cpp +++ b/src/RScore/Individual.cpp @@ -35,11 +35,6 @@ Individual::Individual(Cell* pCell, Patch* pPatch, short stg, short a, short rep float probmale, bool movt, short moveType) { indId = indCounter; indCounter++; // unique identifier for each individual -#if RSDEBUG - //DEBUGLOG << "Individual::Individual(): indId=" << indId - // << " stg=" << stg << " a=" << a << " probmale=" << probmale - // << endl; -#endif fitness = 1.0; stage = stg; if (probmale <= 0.0) sex = FEM; @@ -80,13 +75,6 @@ Individual::Individual(Cell* pCell, Patch* pPatch, short stg, short a, short rep } //pEmigTraits = make_unique(nullptr); //pSettleTraits = make_unique(nullptr); -#if RSDEBUG -//locn currloc = pCurrCell->getLocn(); -//DEBUGLOG << "Individual::Individual(): indId=" << indId -// << " x=" << currloc.x << " y=" << currloc.y -//// << " smsData=" << smsData << " dp=" << smsData->dp -// << endl; -#endif } Individual::~Individual(void) { @@ -281,35 +269,12 @@ void Individual::setSettlementTraits(Species* pSpecies, bool sexDep) { s.beta = getTrait(S_BETA_F)->express(); } -#if RSDEBUG - //DEBUGLOG << "Individual::setSettTraits(): indId=" << indId - // << " s.s0=" << s.s0 << " s.alpha=" << s.alpha << " s.beta=" << s.beta - // << endl; -#endif -#if RSDEBUG - //DEBUGLOG << "Individual::setSettTraits(): indId=" << indId - // << " sparams.s0Mean=" << sparams.s0Mean << " sparams.s0SD=" << sparams.s0SD - // << " sparams.s0Scale=" << sparams.s0Scale - // << endl; -#endif pSettleTraits = make_unique(); pSettleTraits->s0 = (float)(s.s0); pSettleTraits->alpha = (float)(s.alpha); pSettleTraits->beta = (float)(s.beta); -#if RSDEBUG - //DEBUGLOG << "Individual::setSettTraits(): indId=" << indId - // << " setttraits->s0=" << setttraits->s0 - // << " setttraits->alpha=" << setttraits->alpha << " setttraits->beta=" << setttraits->beta - // << endl; -#endif if (pSettleTraits->s0 < 0.0) pSettleTraits->s0 = 0.0; if (pSettleTraits->s0 > 1.0) pSettleTraits->s0 = 1.0; -#if RSDEBUG - //DEBUGLOG << "Individual::setSettTraits(): indId=" << indId - // << " setttraits->s0=" << setttraits->s0 - // << " setttraits->alpha=" << setttraits->alpha << " setttraits->beta=" << setttraits->beta - // << endl; -#endif return; } @@ -371,11 +336,6 @@ void Individual::setYearSteps(int t) { if (t >= 0) path->year = t; else path->year = 666; } -#if RSDEBUG - //DEBUGLOG << "Individual::setYearSteps(): indId=" << indId - // << " t=" << t << " path->year=" << path->year - // << endl; -#endif } pathSteps Individual::getSteps(void) { @@ -429,51 +389,23 @@ void Individual::setEmigTraits(Species* pSpecies, bool sexDep, bool densityDep) } } -#if RSDEBUG - //DEBUGLOG << "Individual::setEmigTraits(): indId=" << indId - // << " eparams.betaMean=" << eparams.betaMean << " eparams.betaSD=" << eparams.betaSD - // << " eparams.betaScale=" << eparams.betaScale - // << endl; -#endif pEmigTraits = make_unique(); pEmigTraits->d0 = (float)(e.d0); pEmigTraits->alpha = (float)(e.alpha); pEmigTraits->beta = (float)(e.beta); -#if RSDEBUG - //DEBUGLOG << "Individual::setEmigTraits(): indId=" << indId - // << " emigtraits->d0=" << emigtraits->d0 - // << " emigtraits->alpha=" << emigtraits->alpha << " emigtraits->beta=" << emigtraits->beta - // << endl; -#endif if (pEmigTraits->d0 < 0.0) pEmigTraits->d0 = 0.0; if (pEmigTraits->d0 > 1.0) pEmigTraits->d0 = 1.0; -#if RSDEBUG - //DEBUGLOG << "Individual::setEmigTraits(): indId=" << indId - // << " emigtraits->d0=" << emigtraits->d0 - // << " emigtraits->alpha=" << emigtraits->alpha << " emigtraits->beta=" << emigtraits->beta - // << endl; -#endif return; } // Get phenotypic emigration traits emigTraits Individual::getEmigTraits(void) { -#if RSDEBUG - //DEBUGLOG << "Individual::getEmigTraits(): indId=" << indId - // << endl; -#endif emigTraits e; e.d0 = e.alpha = e.beta = 0.0; if (pEmigTraits != 0) { e.d0 = pEmigTraits->d0; e.alpha = pEmigTraits->alpha; e.beta = pEmigTraits->beta; } -#if RSDEBUG - //DEBUGLOG << "Individual::getEmigTraits(): indId=" << indId - // << " e.d0=" << e.d0 << " e.alpha=" << e.alpha << " e.beta=" << e.beta - // << endl; -#endif - return e; } // Set phenotypic transfer by kernel traits @@ -502,13 +434,6 @@ void Individual::setKernelTraits(Species* pSpecies, bool sexDep, bool twinKernel float meanDist2 = (float)(k.meanDist2); float probKern1 = (float)(k.probKern1); -#if RSDEBUG - //DEBUGLOG << "Individual::setKernTraits(): indId=" << indId - // << " kerntraits->meanDist1=" << kerntraits->meanDist1 - // << " kerntraits->meanDist2=" << kerntraits->meanDist2 - // << " kerntraits->probKern1=" << kerntraits->probKern1 - // << endl; -#endif if (!pSpecies->useFullKernel()) { // kernel mean(s) may not be less than landscape resolution if (meanDist1 < resol) meanDist1 = (float)resol; @@ -516,13 +441,6 @@ void Individual::setKernelTraits(Species* pSpecies, bool sexDep, bool twinKernel } if (probKern1 < 0.0) probKern1 = 0.0; if (probKern1 > 1.0) probKern1 = 1.0; -#if RSDEBUG - //DEBUGLOG << "Individual::setKernTraits(): indId=" << indId - // << " kerntraits->meanDist1=" << kerntraits->meanDist1 - // << " kerntraits->meanDist2=" << kerntraits->meanDist2 - // << " kerntraits->probKern1=" << kerntraits->probKern1 - // << endl; -#endif auto& pKernel = dynamic_cast(*pTrfrData); pKernel.meanDist1 = meanDist1; pKernel.meanDist2 = meanDist2; @@ -535,10 +453,6 @@ void Individual::setKernelTraits(Species* pSpecies, bool sexDep, bool twinKernel // Get phenotypic emigration traits trfrKernTraits Individual::getKernTraits(void) { -#if RSDEBUG - //DEBUGLOG << "Individual::getKernTraits(): indId=" << indId - // << endl; -#endif trfrKernTraits k; k.meanDist1 = k.meanDist2 = k.probKern1 = 0.0; if (pTrfrData != 0) { @@ -548,12 +462,6 @@ trfrKernTraits Individual::getKernTraits(void) { k.meanDist2 = pKernel.meanDist2; k.probKern1 = pKernel.probKern1; } -#if RSDEBUG - //DEBUGLOG << "Individual::getKernTraits(): indId=" << indId - // << " k.meanDist1=" << k.meanDist1 << " k.meanDist2=" << k.meanDist1 - // << " k.probKern1=" << k.probKern1 - // << endl; -#endif return k; } @@ -571,12 +479,6 @@ void Individual::setSMSTraits(Species* pSpecies) { betaDB = getTrait(SMS_BETADB)->express(); } -#if RSDEBUG - //DEBUGLOG << "Individual::setSMSTraits(): indId=" << indId - // << " dp=" << dp << " gb=" << gb - // << " alphaDB=" << alphaDB << " betaDB=" << betaDB - // << endl; -#endif auto& pSMS = dynamic_cast(*pTrfrData); pSMS.dp = (float)(dp); pSMS.gb = (float)(gb); @@ -588,22 +490,10 @@ void Individual::setSMSTraits(Species* pSpecies) { pSMS.alphaDB = s.alphaDB; pSMS.betaDB = s.betaDB; } -#if RSDEBUG - //DEBUGLOG << "Individual::setSMSTraits() 1111: indId=" << indId - // << " smsData->dp=" << smsData->dp << " smsData->gb=" << smsData->gb - // << " smsData->alphaDB=" << smsData->alphaDB << " smsData->betaDB=" << smsData->betaDB - // << endl; -#endif if (pSMS.dp < 1.0) pSMS.dp = 1.0; if (pSMS.gb < 1.0) pSMS.gb = 1.0; if (pSMS.alphaDB <= 0.0) pSMS.alphaDB = 0.000001f; if (pSMS.betaDB < 1) pSMS.betaDB = 1; -#if RSDEBUG - //DEBUGLOG << "Individual::setSMSTraits() 2222: indId=" << indId - // << " smsData->dp=" << smsData->dp << " smsData->gb=" << smsData->gb - // << " smsData->alphaDB=" << smsData->alphaDB << " smsData->betaDB=" << smsData->betaDB - // << endl; -#endif return; } @@ -613,10 +503,7 @@ trfrData* Individual::getTrfrData(void) { // Get phenotypic transfer by SMS traits trfrSMSTraits Individual::getSMSTraits(void) { -#if RSDEBUG - //DEBUGLOG << "Individual::getSMSTraits(): indId=" << indId << " smsData=" << smsData - // << endl; -#endif + trfrSMSTraits s; s.dp = s.gb = s.alphaDB = 1.0; s.betaDB = 1; if (pTrfrData != 0) { @@ -625,12 +512,7 @@ trfrSMSTraits Individual::getSMSTraits(void) { s.dp = pSMS.dp; s.gb = pSMS.gb; s.alphaDB = pSMS.alphaDB; s.betaDB = pSMS.betaDB; } -#if RSDEBUG - //DEBUGLOG << "Individual::getSMSTraits(): indId=" << indId - // << " s.dp=" << s.dp << " s.gb=" << s.gb - // << " s.alphaDB=" << s.alphaDB << " s.betaDB=" << s.betaDB - // << endl; -#endif + return s; } @@ -647,37 +529,17 @@ void Individual::setCRWTraits(Species* pSpecies, bool sexDep) { c.rho = getTrait(CRW_STEPCORRELATION_F)->express(); } -#if RSDEBUG - //DEBUGLOG << "Individual::setCRWTraits(): indId=" << indId - // << " c.stepLength=" << c.stepLength << " c.rho=" << c.rho - // << endl; -#endif - auto& pCRW = dynamic_cast(*pTrfrData); pCRW.stepLength = (float)(c.stepLength); pCRW.rho = (float)(c.rho); -#if RSDEBUG - //DEBUGLOG << "Individual::setCRWTraits(): indId=" << indId - // << " crw->stepL=" << crw->stepL << " crw->rho=" << crw->rho - // << endl; -#endif if (pCRW.stepLength < 1.0) pCRW.stepLength = 1.0; if (pCRW.rho < 0.0) pCRW.rho = 0.0; if (pCRW.rho > 0.999) pCRW.rho = 0.999f; -#if RSDEBUG - //DEBUGLOG << "Individual::setCRWTraits(): indId=" << indId - // << " crw->stepL=" << crw->stepL << " crw->rho=" << crw->rho - // << endl; -#endif return; } // Get phenotypic transfer by CRW traits trfrCRWTraits Individual::getCRWTraits(void) { -#if RSDEBUG - //DEBUGLOG << "Individual::getCRWTraits(): indId=" << indId - // << endl; -#endif trfrCRWTraits c; c.stepLength = c.rho = 0.0; @@ -688,11 +550,6 @@ trfrCRWTraits Individual::getCRWTraits(void) { c.stepLength = pCRW.stepLength; c.rho = pCRW.rho; } -#if RSDEBUG - //DEBUGLOG << "Individual::getCRWTraits(): indId=" << indId - // << " c.stepLength=" << c.stepLength << " c.rho=" << c.rho - // << endl; -#endif return c; @@ -700,32 +557,16 @@ trfrCRWTraits Individual::getCRWTraits(void) { // Get phenotypic settlement traits settleTraits Individual::getSettTraits(void) { -#if RSDEBUG - //DEBUGLOG << "Individual::getSettTraits(): indId=" << indId - // << endl; -#endif settleTraits s; s.s0 = s.alpha = s.beta = 0.0; if (pSettleTraits != 0) { s.s0 = pSettleTraits->s0; s.alpha = pSettleTraits->alpha; s.beta = pSettleTraits->beta; } -#if RSDEBUG - //DEBUGLOG << "Individual::getSettTraits(): indId=" << indId - // << " s.s0=" << s.s0 << " s.alpha=" << s.alpha << " s.beta=" << s.beta - // << endl; -#endif return s; } -/* -locus Individual::getAlleles(int g) { -locus l; l.allele[0] = l.allele[1] = 0.0; -if (pGenome != 0) l = pGenome->getAlleles(g); -return l; -} -*/ void Individual::setStatus(short s) { if (s >= 0 && s <= 9) status = s; @@ -827,19 +668,8 @@ int Individual::moveKernel(Landscape* pLandscape, Species* pSpecies, } } } -#if RSDEBUG - //Patch *startPatch = (Patch*)startpatch; - //DEBUGLOG << "Individual::moveKernel(): indId=" << indId << " x=" << loc.x << " y=" << loc.y - //// << " natalPatch = " << natalPatch - //// << " startpatch = " << startpatch << " patchNum = " << startPatch->getPatchNum() - // << " kern.meanDist1=" << kern.meanDist1; - //if (trfr.twinKern) { - // DEBUGLOG << " meanDist2=" << kern.meanDist2 << " probKern1=" << kern.probKern1; - //} - //DEBUGLOG << endl; -#endif -// scale the appropriate kernel mean to the cell size + // scale the appropriate kernel mean to the cell size if (trfr.twinKern) { if (pRandom->Bernoulli(kern.probKern1)) @@ -849,27 +679,9 @@ int Individual::moveKernel(Landscape* pLandscape, Species* pSpecies, } else meandist = kern.meanDist1 / (float)land.resol; -#if RSDEBUG - //DEBUGLOG << "Individual::moveKernel(): indId=" << indId << " meandist=" << meandist << endl; -#endif -// scaled mean may not be less than 1 unless emigration derives from the kernel -// (i.e. the 'use full kernel' option is applied) + // scaled mean may not be less than 1 unless emigration derives from the kernel + // (i.e. the 'use full kernel' option is applied) if (!usefullkernel && meandist < 1.0) meandist = 1.0; -#if RSDEBUG - //DEBUGLOG << "Individual::moveKernel(): indId=" << indId << " meandist=" << meandist << endl; -#endif - -#if RSDEBUG -//Patch *startPatch = (Patch*)startpatch; -//DEBUGLOG << "Individual::moveKernel(): indId = " << indId << " x = " << x << " y = " << y -// << " natalPatch = " << natalPatch -//// << " startpatch = " << startpatch << " patchNum = " << startPatch->getPatchNum() -// << " meanDist1 = " << kern.meanDist1; -//if (trfr.twinKern) { -// DEBUGLOG << " probKern1 = " << kern.probKern1 << " meanDist2 = " << kern.meanDist2; -//} -//DEBUGLOG << " meandist = " << meandist << endl; -#endif int loopsteps = 0; // new counter to prevent infinite loop added 14/8/15 do { @@ -904,12 +716,6 @@ int Individual::moveKernel(Landscape* pLandscape, Species* pSpecies, if (path != 0) (path->year)++; #endif loopsteps++; -#if RSDEBUG - //DEBUGLOG << "Individual::moveKernel(): indId=" << indId << " status=" << status - // << " loopsteps=" << loopsteps << " newX=" << newX << " newY=" << newY - // << " loc.x=" << loc.x << " loc.y=" << loc.y - // << endl; -#endif } while (loopsteps < 1000 && ((!absorbing && (newX < land.minX || newX > land.maxX || newY < land.minY || newY > land.maxY)) @@ -945,12 +751,6 @@ int Individual::moveKernel(Landscape* pLandscape, Species* pSpecies, patch = 0; patchNum = -1; } -#if RSDEBUG - //DEBUGLOG << "Individual::moveKernel(): indId=" << indId << " status=" << status - // << " loopsteps=" << loopsteps << " newX=" << newX << " newY=" << newY - // << " pCell=" << pCell << " patch=" << patch << " patchNum=" << patchNum - // << endl; -#endif } while (!absorbing && patchNum < 0 && loopsteps < 1000); // in a no-data region } while (!usefullkernel && pPatch == pNatalPatch && loopsteps < 1000); // still in the original (natal) patch @@ -987,15 +787,8 @@ int Individual::moveKernel(Landscape* pLandscape, Species* pSpecies, status = 6; dispersing = 0; } -#if RSDEBUG - //DEBUGLOG << "Individual::moveKernel(): indId=" << indId - // << " newX=" << newX << " newY=" << newY - // << " patch=" << patch - // << " patchNum=" << patchNum << " status=" << status; - //DEBUGLOG << endl; -#endif -// apply dispersal-related mortality, which may be distance-dependent + // apply dispersal-related mortality, which may be distance-dependent dist *= (float)land.resol; // re-scale distance moved to landscape scale if (status < 7) { double dispmort; @@ -1044,11 +837,6 @@ int Individual::moveStep(Landscape* pLandscape, Species* pSpecies, settleSteps settsteps = pSpecies->getSteps(stage, sex); patch = pCurrCell->getPatch(); -#if RSDEBUG - //DEBUGLOG << "Individual::moveStep() AAAA: indId=" << indId - // << " pCurrCell=" << pCurrCell << " patch=" << patch - // << endl; -#endif if (patch == 0) { // matrix pPatch = 0; @@ -1066,35 +854,12 @@ int Individual::moveStep(Landscape* pLandscape, Species* pSpecies, mortprob = 1.0; } else mortprob = pSpecies->getHabMort(h); -#if RSDEBUG - //locn temploc = pCurrCell->getLocn(); - //DEBUGLOG << "Individual::moveStep(): x=" << temploc.x << " y=" << temploc.x - // << " landIx=" << landIx << " h=" << h << " mortprob=" << mortprob - // << endl; -#endif } else mortprob = movt.stepMort; // ... unless individual has not yet left natal patch in emigration year if (pPatch == pNatalPatch && path->out == 0 && path->year == path->total) { mortprob = 0.0; } -#if RSDEBUG - locn loc0, loc1, loc2; - //loc0 = pCurrCell->getLocn(); - //DEBUGLOG << "Individual::moveStep() BBBB: indId=" << indId << " status=" << status - // << " path->year=" << path->year << " path->out=" << path->out - // << " settleStatus=" << path->settleStatus - // << " x=" << loc0.x << " y=" << loc0.y - //// << " patch=" << patch - // << " pPatch=" << pPatch - // << " patchNum=" << patchNum; - //// << " natalPatch=" << natalPatch; - ////if (crw != 0) { - //// DEBUGLOG << " xc=" << crw->xc << " yc=" << crw->yc; - //// DEBUGLOG << " rho=" << movt.rho << " stepLength=" << movt.stepLength; - ////} - //DEBUGLOG << endl; -#endif if (pRandom->Bernoulli(mortprob)) { // individual dies status = 7; dispersing = 0; @@ -1113,23 +878,7 @@ int Individual::moveStep(Landscape* pLandscape, Species* pSpecies, switch (trfr.moveType) { case 1: // SMS -#if RSDEBUG - //loc1 = pCurrCell->getLocn(); - //DEBUGLOG << "Individual::moveStep() FFFF: indId=" << indId << " status=" << status - //// << " path->year=" << path->year - // << " path->season=" << path->season - // << " x=" << loc1.x << " y=" << loc1.y - // << " smsData->goalType=" << smsData->goalType - // << " goal.x=" << smsData->goal.x - // << " goal.y=" << smsData->goal.y - // << endl; -#endif move = smsMove(pLandscape, pSpecies, landIx, pPatch == pNatalPatch, trfr.indVar, absorbing); -#if RSDEBUG - //DEBUGLOG << "Individual::moveStep() GGGG: indId=" << indId << " status=" << status - // << " move.dist=" << move.dist - // << endl; -#endif if (move.dist < 0.0) { // either INTERNAL ERROR CONDITION - INDIVIDUAL IS IN NO-DATA SQUARE // or individual has crossed absorbing boundary ... @@ -1138,16 +887,8 @@ int Individual::moveStep(Landscape* pLandscape, Species* pSpecies, dispersing = 0; } else { -#if RSDEBUG - //loc1 = pCurrCell->getLocn(); - //DEBUGLOG << "Individual::moveStep() HHHH: indId=" << indId << " status=" << status - // << " path->year=" << path->year - // << " x=" << loc1.x << " y=" << loc1.y - //// << " smsData = " << smsData - // << endl; -#endif - // WOULD IT BE MORE EFFICIENT FOR smsMove TO RETURN A POINTER TO THE NEW CELL? ... + // WOULD IT BE MORE EFFICIENT FOR smsMove TO RETURN A POINTER TO THE NEW CELL? ... patch = pCurrCell->getPatch(); //int patchnum; @@ -1169,7 +910,7 @@ int Individual::moveStep(Landscape* pLandscape, Species* pSpecies, auto & pCRW = dynamic_cast(*pTrfrData); - if (trfr.indVar) { + if (trfr.indVar) { movt.stepLength = pCRW.stepLength; movt.rho = pCRW.rho; } @@ -1203,13 +944,6 @@ int Individual::moveStep(Landscape* pLandscape, Species* pSpecies, if (xcnew < 0.0) newX = -1; else newX = (int)xcnew; if (ycnew < 0.0) newY = -1; else newY = (int)ycnew; loopsteps++; -#if RSDEBUG - //DEBUGLOG << "Individual::moveStep(): indId=" << indId - // << " xc=" << crw->xc << " yc=" << crw->yc << " pCurrCell=" << pCurrCell - // << " steps=" << path->year << " loopsteps=" << loopsteps - // << " steplen=" << steplen << " rho=" << rho << " angle=" << angle - // << " xcnew=" << xcnew << " ycnew=" << ycnew << " newX=" << newX << " newY=" << newY << endl; -#endif } while (!absorbing && loopsteps < 1000 && (newX < land.minX || newX > land.maxX || newY < land.minY || newY > land.maxY)); if (newX < land.minX || newX > land.maxX || newY < land.minY || newY > land.maxY) @@ -1222,11 +956,6 @@ int Individual::moveStep(Landscape* pLandscape, Species* pSpecies, } else patch = pCurrCell->getPatch(); -#if RSDEBUG - //DEBUGLOG << "Individual::moveStep(): indId=" << indId - // << " loopsteps=" << loopsteps << " absorbed=" << absorbed - // << " pCurrCell=" << pCurrCell << " patch=" << patch << endl; -#endif } while (!absorbing && pCurrCell == 0 && loopsteps < 1000); pCRW.prevdrn = (float)angle; pCRW.xc = (float)xcnew; pCRW.yc = (float)ycnew; @@ -1246,36 +975,10 @@ int Individual::moveStep(Landscape* pLandscape, Species* pSpecies, pCurrCell = pPrevCell; } } -#if RSDEBUG - //DEBUGLOG << "Individual::moveStep(): indId=" << indId - // << " status=" << status - // << " pCurrCell=" << pCurrCell << " patch=" << patch << endl; -#endif break; } // end of switch (trfr.moveType) -#if RSDEBUG -//locn loc2; -//if (pCurrCell > 0) { -// loc2 = pCurrCell->getLocn(); -//} -//else { -// loc2.x = -9999; loc2.y = -9999; -//} -//DEBUGLOG << "Individual::moveStep() ZZZZ: indId=" << indId -// << " status=" << status -// << " path->total=" << path->total -// << " x=" << loc2.x << " y=" << loc2.y -// << " patch=" << patch; -//if (patch > 0) { -// pPatch = (Patch*)patch; -// DEBUGLOG << " patchNum=" << pPatch->getPatchNum() -// << " getK()=" << pPatch->getK() -// << " popn=" << pPatch->getPopn((int)pSpecies); -//} -// DEBUGLOG << endl; -#endif if (patch > 0 // not no-data area or matrix && path->total >= settsteps.minSteps) { pPatch = (Patch*)patch; @@ -1324,12 +1027,6 @@ movedata Individual::smsMove(Landscape* pLand, Species* pSpecies, locn current; auto& pSMS = dynamic_cast(*pTrfrData); - //if (write_out) { - // out<findCell(x,y); if (pCurrCell == 0) { // x,y is a NODATA square - this should not occur here @@ -1338,16 +1035,11 @@ movedata Individual::smsMove(Landscape* pLand, Species* pSpecies, return move; } -#if RSDEBUG - //DEBUGLOG << "Individual::smsMove(): this=" << this << endl; -#endif - landData land = pLand->getLandData(); trfrSMSTraits movt = pSpecies->getSMSTraits(); current = pCurrCell->getLocn(); //get weights for directional persistence.... - //if ((path->out > 0 && path->out < 10 && path->out < 2*movt.pr) if ((path->out > 0 && path->out <= (movt.pr + 1)) || natalPatch || (movt.straigtenPath && path->settleStatus > 0)) { @@ -1360,25 +1052,8 @@ movedata Individual::smsMove(Landscape* pLand, Species* pSpecies, else nbr = getSimDir(current.x, current.y, movt.dp); } if (natalPatch || path->settleStatus > 0) path->out = 0; - //if (natalPatch) path->out = 0; -#if RSDEBUG -//DEBUGLOG << "Individual::smsMove() 0000: nbr matrix" << endl; -//for (y2 = 2; y2 > -1; y2--) { -// for (x2 = 0; x2 < 3; x2++) DEBUGLOG << nbr.cell[x2][y2] << " "; -// DEBUGLOG << endl; -//} -#endif -//if (write_out) { -// out< -1; y2--) { -// for (x2 = 0; x2 < 3; x2++) out< -1; y2--) { - // for (x2 = 0; x2 < 3; x2++) out< -1; y2--) { -// for (x2 = 0; x2 < 3; x2++) out< -1; y2--) { - // for (x2 = 0; x2 < 3; x2++) DEBUGLOG << hab.cell[x2][y2] << " "; - // DEBUGLOG << endl; - //} -#endif pCurrCell->setEffCosts(hab); } else { // they have already been calculated - no action required - // if (write_out) { - // out<<"*** using previous effective costs ***"< -1; y2--) { - // for (x2 = 0; x2 < 3; x2++) { - // out< -1; y2--) { - // for (x2 = 0; x2 < 3; x2++) DEBUGLOG << nbr.cell[x2][y2] << " "; - // DEBUGLOG << endl; - //} -#endif -// determine reciprocal of effective cost for the 8 neighbours -//if (write_out) out<<"reciprocal weighted effective costs:"< -1; y2--) { for (x2 = 0; x2 < 3; x2++) { if (nbr.cell[x2][y2] > 0.0) nbr.cell[x2][y2] = 1.0f / nbr.cell[x2][y2]; - // if (write_out) { - // out< -1; y2--) { -// for (x2 = 0; x2 < 3; x2++) { -// temp.cell[x2][y2] = nbr.cell[x2][y2]; -// if (current.x == 488 && current.y == 422) { -// pCell = pLand->findCell((current.x+x2-1),(current.y+y2-1)); -// DEBUGLOG << "Individual::smsMove(): this=" << this -// << " IN THE PROBLEM CELL" -// << " y=" << current.y << " x=" << current.x -// << " y2=" << y2 << " x2=" << x2 -// << " pCell=" << pCell; -// if (pCell != 0) DEBUGLOG << " pCell->getCost=" << pCell->getCost(); -// DEBUGLOG << endl; -// } -// } -//} -#endif - for (y2 = 2; y2 > -1; y2--) { for (x2 = 0; x2 < 3; x2++) { if (!absorbing) { @@ -1537,44 +1128,20 @@ movedata Individual::smsMove(Landscape* pLand, Species* pSpecies, if (pCell == 0) nbr.cell[x2][y2] = 0.0; // no-data cell } } -#if RSDEBUG - //DEBUGLOG << "Individual::smsMove(): this=" << this - // << " y=" << current.y << " x=" << current.x - // << " y2=" << y2 << " x2=" << x2 - // << " pCell=" << pCell - // << endl; -#endif -// if (write_out) { -// out< 0.0) { // should always be the case, but safest to check... for (y2 = 2; y2 > -1; y2--) { for (x2 = 0; x2 < 3; x2++) { nbr.cell[x2][y2] = nbr.cell[x2][y2] / (float)sum_nbrs; - // if (write_out) { - // out< -1; y2--) { - // for (x2 = 0; x2 < 3; x2++) DEBUGLOG << nbr.cell[x2][y2] << " "; - // DEBUGLOG << endl; - //} -#endif -// set up cell selection probabilities -//if (write_out) out<<"rnd = "<findCell(newX, newY); } } while (!absorbing && pNewCell == 0 && loopsteps < 1000); // no-data cell -#if RSDEBUG - //DEBUGLOG << "Individual::smsMove() 8888: pNewCell=" << pNewCell - // << " loopsteps=" << loopsteps - // << " current.x=" << current.x << " current.y=" << current.y - // << " newX=" << newX << " newY=" << newY - // << " land.minX=" << land.minX << " land.minY=" << land.minY - // << " land.maxX=" << land.maxX << " land.maxY=" << land.maxY - // << endl; -#endif if (loopsteps >= 1000 || pNewCell == 0) { // unable to make a move or crossed absorbing boundary // flag individual to die @@ -1688,15 +1233,12 @@ array3x3d Individual::getSimDir(const int x, const int y, const float dp) if ((x - prev.x) == 0 && (y - prev.y) == 0) { // back to 'square 1' (first memory location) - use previous step drn only prev = memory.back(); - // if (write_out) out<<"step 3"<getHabIndex(landIx); cost = pSpecies->getHabCost(h); -#if RSDEBUG - //DEBUGLOG << "Individual::getHabMatrix(): x4=" << x4 << " y4=" << y4 - // << " landIx=" << landIx << " h=" << h << " cost=" << cost - // << endl; -#endif pCell->setCost(cost); } else { -#if RSDEBUG - //DEBUGLOG << "Individual::getHabMatrix(): x4=" << x4 << " y4=" << y4 - // << " cost=" << cost - // << endl; -#endif - + // nothing? } } } @@ -1938,26 +1457,14 @@ array3x3f Individual::getHabMatrix(Landscape* pLand, Species* pSpecies, ncells++; sumweights += weight; } } - // if (write_out2) out2<get_target() > 50) targetseen++; - // } - //#endif } //end of y3 loop } //end of x3 loop - // if (write_out) out<<"ncells in PR = "<(copyFrom); + const crwData& pCopy = dynamic_cast(copyFrom); stepLength = pCopy.stepLength; rho = pCopy.rho; @@ -155,7 +154,7 @@ struct smsData : trfrData { //static float stepMort; //static bool straigtenPath; - smsData(locn prevA, locn goalA) : prev(prevA), goal(goalA), dp(0.0), gb(0.0), alphaDB(0.0), betaDB(0.0) {} + smsData(locn prevA, locn goalA) : prev(prevA), goal(goalA), dp(0.0), gb(0.0), alphaDB(0.0), betaDB(0) {} ~smsData() {} @@ -207,7 +206,7 @@ struct kernelData : trfrData { movement_t getType() { return KERNEL; } void clone(const trfrData& copyFrom) { - auto pCopy = dynamic_cast(copyFrom); + const kernelData& pCopy = dynamic_cast(copyFrom); meanDist1 = pCopy.meanDist1; meanDist2 = pCopy.meanDist2; probKern1 = pCopy.probKern1; @@ -334,13 +333,6 @@ class Individual { const short, // landscape change index const bool // absorbing boundaries? ); - void drawMove( // Visualise paths resulting from movement simulation model - // NULL for the batch version - const float, // initial x co-ordinate - const float, // initial y co-ordinate - const float, // final x co-ordinate - const float // final y co-ordinate - ); movedata smsMove( // Move to a neighbouring cell according to the SMS algorithm Landscape*, // pointer to Landscape Species*, // pointer to Species @@ -435,5 +427,9 @@ extern ofstream DEBUGLOG; extern ofstream outMovePaths; #endif -//--------------------------------------------------------------------------- +#if RSDEBUG +void testIndividual(); #endif + +//--------------------------------------------------------------------------- +#endif // IndividualH diff --git a/src/RScore/Landscape.cpp b/src/RScore/Landscape.cpp index e38d93e..369c3a3 100644 --- a/src/RScore/Landscape.cpp +++ b/src/RScore/Landscape.cpp @@ -208,10 +208,6 @@ int InitDist::readDistribution(string distfile) { if (dfile >> p) { #else dfile >> p; -#endif -#if RSDEBUG - //DEBUGLOG << "InitDist::readDistribution():" - // << " y = " << y << " x = " << x << " p = " << p << endl; #endif if (p == nodata || p == 0 || p == 1) { // only valid values if (p == 1) { // species present @@ -266,60 +262,26 @@ Landscape::Landscape(void) { pix = (int)gpix; minEast = minNorth = 0.0; cells = 0; -#if RSDEBUG - // NOTE: do NOT write to output stream before it has been opened - it will be empty - //DebugGUI("Landscape::Landscape(): this = " + Int2Str((int)this) - // + " pCell = " + Int2Str((int)pCell)); - //DEBUGLOGGUI << "Landscape::Landscape(): this = " << this << " pCell = " << pCell << endl; -#endif connectMatrix = 0; epsGlobal = 0; patchChgMatrix = 0; costsChgMatrix = 0; - -#if RSDEBUG - //DEBUGLOG << "Landscape::Landscape():" - // << " rasterType= " << rasterType << endl; -#endif -#if RSDEBUG -//MemoLine(("Landscape::Landscape(): landscape created " + Int2Str(0) -// ).c_str()); -#endif } Landscape::~Landscape() { -#if RSDEBUG - //DebugGUI(("Landscape::~Landscape(): this=" + Int2Str((int)this) + " cells=" + Int2Str((int)cells) - // + " maxX=" + Int2Str(maxX) + " maxY=" + Int2Str(maxY)).c_str()); -#endif if (cells != 0) { for (int y = dimY - 1; y >= 0; y--) { -#if RSDEBUG - //DebugGUI(("Landscape::~Landscape(): y=" + Int2Str(y) + " cells[y]=" + Int2Str((int)cells[y])).c_str()); -#endif for (int x = 0; x < dimX; x++) { -#if RSDEBUG - //DebugGUI(("Landscape::~Landscape(): y=" + Int2Str(y) + " x=" + Int2Str(x) - // + " cells[y][x]=" + Int2Str((int)cells[y][x])).c_str()); -#endif if (cells[y][x] != 0) delete cells[y][x]; } if (cells[y] != 0) { -#if RSDEBUG - //DebugGUI(("Landscape::~Landscape(): deleting cells[y]=" + Int2Str((int)cells[y])).c_str()); -#endif -// delete cells[y]; delete[] cells[y]; } } delete[] cells; cells = 0; } -#if RSDEBUG - //DebugGUI(("Landscape::~Landscape(): this=" + Int2Str((int)this) - // + " cells=" + Int2Str((int)cells)).c_str()); -#endif int npatches = (int)patches.size(); for (int i = 0; i < npatches; i++) @@ -345,63 +307,28 @@ Landscape::~Landscape() { deleteConnectMatrix(); deletePatchChgMatrix(); if (epsGlobal != 0) delete[] epsGlobal; - -#if RSDEBUG - //MemoLine(("Landscape::~Landscape(): landscape deleted " + Int2Str(1) - // ).c_str()); -#endif } // Remove all patches and cells // Used for replicating artificial landscape without deleting the landscape itself void Landscape::resetLand(void) { - -#if RSDEBUG - //DebugGUI("Landscape::resetLand(): starting..."); -#endif resetLandLimits(); int npatches = (int)patches.size(); -#if RSDEBUG - //DebugGUI(("Landscape::resetLand(): npatches=" + Int2Str(npatches) - // ).c_str()); -#endif for (int i = 0; i < npatches; i++) if (patches[i] != NULL) delete patches[i]; patches.clear(); patchnums.clear(); -#if RSDEBUG - //DebugGUI("Landscape::resetLand(): finished resetting patches"); -#endif - -#if RSDEBUG -//DebugGUI(("Landscape::resetLand(): this=" + Int2Str((int)this) -// + " cells=" + Int2Str((int)cells)).c_str()); -#endif if (cells != 0) { for (int y = dimY - 1; y >= 0; y--) { -#if RSDEBUG - //DebugGUI(("Landscape::resetLand(): y=" + Int2Str(y) + " cells[y]=" + Int2Str((int)cells[y])).c_str()); -#endif for (int x = 0; x < dimX; x++) { -#if RSDEBUG - //DebugGUI(("Landscape::resetLand(): y=" + Int2Str(y) + " x=" + Int2Str(x) - // + " cells[y][x]=" + Int2Str((int)cells[y][x])).c_str()); -#endif if (cells[y][x] != 0) delete cells[y][x]; } if (cells[y] != 0) { -#if RSDEBUG - //DebugGUI(("Landscape::resetLand(): deleting cells[y]=" + Int2Str((int)cells[y])).c_str()); -#endif delete[] cells[y]; } } delete[] cells; cells = 0; } -#if RSDEBUG - //DebugGUI(("Landscape::resetLand(): this=" + Int2Str((int)this) - // + " cells=" + Int2Str((int)cells)).c_str()); -#endif } void Landscape::setLandParams(landParams ppp, bool batchMode) @@ -593,30 +520,15 @@ int Landscape::colourCount(void) { //--------------------------------------------------------------------------- void Landscape::setCellArray(void) { -#if RSDEBUG - //DebugGUI(("Landscape::setCellArray(): start: this=" + Int2Str((int)this) - // + " cells=" + Int2Str((int)cells)).c_str()); -#endif if (cells != 0) resetLand(); //cells = new Cell **[maxY+1]; cells = new Cell * *[dimY]; -#if RSDEBUG - //DebugGUI(("Landscape::setCellArray(): cells=" + Int2Str((int)cells)).c_str()); -#endif for (int y = dimY - 1; y >= 0; y--) { cells[y] = new Cell * [dimX]; -#if RSDEBUG - //DebugGUI(("Landscape::setCellArray(): y=" + Int2Str(y) - // + " cells[y]=" + Int2Str((int)cells[y])).c_str()); -#endif for (int x = 0; x < dimX; x++) { cells[y][x] = 0; } } -#if RSDEBUG - //DebugGUI(("Landscape::setCellArray(): end: this=" + Int2Str((int)this) - // + " cells=" + Int2Str((int)cells)).c_str()); -#endif } void Landscape::addPatchNum(int p) { @@ -635,7 +547,7 @@ void Landscape::addPatchNum(int p) { /* Create an artificial landscape (random or fractal), which can be either binary (habitat index 0 is the matrix, 1 is suitable habitat) or continuous (0 is the matrix, >0 is suitable habitat) */ -void Landscape::generatePatches(void) +void Landscape::generatePatches(Species* pSpecies) { int x, y, ncells; double p; @@ -644,17 +556,6 @@ void Landscape::generatePatches(void) vector ArtLandscape; -#if RSDEBUG - //int iiiiii = (int)fractal; - //DEBUGLOG << "Landscape::generatePatches(): (int)fractal=" << iiiiii - // << " rasterType=" << rasterType - // << " continuous=" << continuous - // << " dimX=" << dimX << " dimY=" << dimY - // << " propSuit=" << propSuit - // << " hurst=" << hurst - // << " maxPct=" << maxPct << " minPct=" << minPct - // << endl; -#endif setCellArray(); @@ -668,9 +569,6 @@ void Landscape::generatePatches(void) // habitat cells for (int yy = dimY - 1; yy >= 0; yy--) { for (int xx = 0; xx < dimX; xx++) { -#if RSDEBUG - //DEBUGLOG << "Landscape::generatePatches(): yy=" << yy << " xx=" << xx << endl; -#endif addNewCellToLand(xx, yy, 0); } } @@ -689,11 +587,6 @@ void Landscape::generatePatches(void) vector::iterator iter = ArtLandscape.begin(); while (iter != ArtLandscape.end()) { x = iter->y_coord; y = iter->x_coord; -#if RSDEBUG - //DEBUGLOG << "Landscape::generatePatches(): x=" << x << " y=" << y - // << " iter->avail=" << iter->avail - // << " iter->value=" << iter->value << endl; -#endif pCell = findCell(x, y); if (continuous) { if (iter->value > 0.0) { // habitat @@ -722,12 +615,6 @@ void Landscape::generatePatches(void) else { // random landscape int hab = 0; ncells = (int)((float)(dimX) * (float)(dimY)*propSuit + 0.00001); // no. of cells to initialise -#if RSDEBUG - //DEBUGLOG << "Landscape::generatePatches(): dimX=" << dimX << " dimY=" << dimY - // << " propSuit=" << propSuit - // << " PRODUCT=" << ((float)(dimX) * (float)(dimY) * propSuit + 0.00001) - // << " ncells=" << ncells << endl; -#endif int i = 0; do { do { @@ -735,11 +622,6 @@ void Landscape::generatePatches(void) pCell = findCell(x, y); hab = pCell->getHabIndex(0); } while (hab > 0); -#if RSDEBUG - //DEBUGLOG << "Landscape::generatePatches() 00000: y=" << y << " x=" << x - // << " i=" << i << " hab=" << hab << " patchnum=" << patchnum - // << endl; -#endif patchnums.push_back(patchnum); pPatch = newPatch(patchnum++); pCell = findCell(x, y); @@ -764,10 +646,6 @@ void Landscape::generatePatches(void) } else { // discrete if (pCell->getHabIndex(0) == 0) { -#if RSDEBUG - //DEBUGLOG << "Landscape::generatePatches() 11111: yy=" << yy << " xx=" << xx - // << endl; -#endif addCellToPatch(pCell, patches[0], x); } } @@ -775,11 +653,6 @@ void Landscape::generatePatches(void) } } -#if RSDEBUG - //DEBUGLOG << "Landscape::generatePatches(): finished, no. of patches = " - // << patchCount() << endl; -#endif - simParams sim = paramsSim->getSim(); if (pSpecies->getNumberOfNeutralLoci() > 0 && (sim.outputWCFstat || sim.outputPairwiseFst)) { @@ -802,11 +675,6 @@ void Landscape::generatePatches(void) habitat cells are added to the matrix patch) */ void Landscape::allocatePatches(Species* pSpecies) { -#if RSDEBUG - //DEBUGLOG << "Landscape::allocatePatches(): pSpecies=" << pSpecies - // << " N patches=" << (int)patches.size() << endl; -#endif -//int hx; float habK; Patch* pPatch; Cell* pCell; @@ -814,10 +682,6 @@ void Landscape::allocatePatches(Species* pSpecies) // delete all existing patches int npatches = (int)patches.size(); for (int i = 0; i < npatches; i++) { -#if RSDEBUG - //DEBUGLOG << "Landscape::allocatePatches(): i=" << i - // << " patches[i]=" << patches[i] << endl; -#endif if (patches[i] != NULL) delete patches[i]; } patches.clear(); @@ -825,10 +689,6 @@ void Landscape::allocatePatches(Species* pSpecies) patches.push_back(new Patch(0, 0)); Patch* matrixPatch = patches[0]; patchnums.push_back(0); -#if RSDEBUG - //DEBUGLOG << "Landscape::allocatePatches(): npatches=" << npatches - // << " N patches=" << (int)patches.size() << endl; -#endif int patchnum = 1; switch (rasterType) { @@ -836,10 +696,6 @@ void Landscape::allocatePatches(Species* pSpecies) case 0: // habitat codes for (int y = dimY - 1; y >= 0; y--) { for (int x = 0; x < dimX; x++) { -#if RSDEBUG - //DEBUGLOG << "Landscape::allocatePatches(): x=" << x << " y=" << y - // << " cells=" << cells[y][x] << endl; -#endif if (cells[y][x] != 0) { // not no-data cell pCell = cells[y][x]; // hx = pCell->getHabIndex(); @@ -848,13 +704,6 @@ void Landscape::allocatePatches(Species* pSpecies) int nhab = pCell->nHabitats(); for (int i = 0; i < nhab; i++) { habK += pSpecies->getHabK(pCell->getHabIndex(i)); -#if RSDEBUG - //DEBUGLOG << "Landscape::allocatePatches(): x=" << x << " y=" << y - // << " i=" << i - // << " cells[y][x]=" << cells[y][x] - // << " pCell=" << pCell - // << " habK=" << habK << endl; -#endif } if (habK > 0.0) { // cell is suitable - create a patch for it patchnums.push_back(patchnum); @@ -872,10 +721,6 @@ void Landscape::allocatePatches(Species* pSpecies) case 1: // habitat cover for (int y = dimY - 1; y >= 0; y--) { for (int x = 0; x < dimX; x++) { -#if RSDEBUG - //DEBUGLOG << "Landscape::allocatePatches(): x=" << x << " y=" << y - // << " cells=" << cells[y][x] << endl; -#endif if (cells[y][x] != 0) { // not no-data cell pCell = cells[y][x]; habK = 0.0; @@ -884,13 +729,6 @@ void Landscape::allocatePatches(Species* pSpecies) for (int i = 0; i < nhab; i++) { habK += pSpecies->getHabK(i) * pCell->getHabitat(i) / 100.0f; -#if RSDEBUG - //DEBUGLOG << "Landscape::allocatePatches(): x=" << x << " y=" << y - // << " i=" << i - // << " cells[y][x]=" << cells[y][x] - // << " pCell=" << pCell - // << " habK=" << habK << endl; -#endif } if (habK > 0.0) { // cell is suitable - create a patch for it patchnums.push_back(patchnum); @@ -908,25 +746,13 @@ void Landscape::allocatePatches(Species* pSpecies) case 2: // habitat quality for (int y = dimY - 1; y >= 0; y--) { for (int x = 0; x < dimX; x++) { -#if RSDEBUG - //DEBUGLOG << "Landscape::allocatePatches(): x=" << x << " y=" << y - // << " cells=" << cells[y][x] << endl; -#endif if (cells[y][x] != 0) { // not no-data cell pCell = cells[y][x]; habK = 0.0; int nhab = pCell->nHabitats(); - // for (int i = 0; i < nHab; i++) for (int i = 0; i < nhab; i++) { habK += pSpecies->getHabK(0) * pCell->getHabitat(i) / 100.0f; -#if RSDEBUG - //DEBUGLOG << "Landscape::allocatePatches(): x=" << x << " y=" << y - // << " i=" << i - // << " cells[y][x]=" << cells[y][x] - // << " pCell=" << pCell - // << " habK=" << habK << endl; -#endif } if (habK > 0.0) { // cell is suitable (at some time) - create a patch for it patchnums.push_back(patchnum); @@ -955,12 +781,6 @@ void Landscape::allocatePatches(Species* pSpecies) pSpecies->setSamplePatchList(patchesToSample); } - - -#if RSDEBUG - //DEBUGLOG << "Landscape::allocatePatches(): finished, N patches = " << (int)patches.size() << endl; -#endif - } Patch* Landscape::newPatch(int num) @@ -1118,12 +938,6 @@ void Landscape::updateCarryingCapacity(Species* pSpecies, int yr, short landIx) patchLimits landlimits; landlimits.xMin = minX; landlimits.xMax = maxX; landlimits.yMin = minY; landlimits.yMax = maxY; -#if RSDEBUG - //DEBUGLOG << "Landscape::updateCarryingCapacity(): yr=" << yr - // << " xMin=" << landlimits.xMin << " yMin=" << landlimits.yMin - // << " xMax=" << landlimits.xMax << " yMax=" << landlimits.yMax - // << endl; -#endif int npatches = (int)patches.size(); for (int i = 0; i < npatches; i++) { if (patches[i]->getPatchNum() != 0) { // not matrix patch @@ -1181,13 +995,6 @@ int Landscape::checkTotalCover(void) { float sumCover = 0.0; for (int i = 0; i < nHab; i++) { sumCover += cells[y][x]->getHabitat(i); -#if RSDEBUG - //DebugGUI(("Landscape::checkTotalCover(): y=" + Int2Str(y) - // + " x=" + Int2Str(x) - // + " i=" + Int2Str(i) - // + " sumCover=" + Float2Str(sumCover) - // ).c_str()); -#endif } if (sumCover > 100.00001) nCells++; // decimal part to allow for floating point error if (sumCover <= 0.0) // cell is a matrix cell @@ -1209,28 +1016,13 @@ void Landscape::updateHabitatIndices(void) { // convert codes in landscape int h; int changes = (int)landchanges.size(); -#if RSDEBUG - //DebugGUI(("Landscape::updateHabitatIndices(): nHab=" + Int2Str(nHab) - // + " changes=" + Int2Str(changes) - // ).c_str()); -#endif for (int y = dimY - 1; y >= 0; y--) { for (int x = 0; x < dimX; x++) { if (cells[y][x] != 0) { // not a no-data cell for (int c = 0; c <= changes; c++) { h = cells[y][x]->getHabIndex(c); -#if RSDEBUG - //DebugGUI(("Landscape::updateHabitatIndices() 00000: y=" + Int2Str(y) - // + " x=" + Int2Str(x) + " c=" + Int2Str(c) + " h=" + Int2Str(h) - // ).c_str()); -#endif if (h >= 0) { h = findHabCode(h); -#if RSDEBUG - //DebugGUI(("Landscape::updateHabitatIndices() 11111: y=" + Int2Str(y) - // + " x=" + Int2Str(x) + " c=" + Int2Str(c) + " h=" + Int2Str(h) - // ).c_str()); -#endif cells[y][x]->changeHabIndex(c, h); } } @@ -1248,11 +1040,6 @@ void Landscape::setEnvGradient(Species* pSpecies, bool initial) double envval; // gradient parameters envGradParams grad = paramsGrad->getGradient(); -#if RSDEBUG - //DEBUGLOG << "Landscape::setEnvGradient(): grad.opt_y = " << grad.opt_y - // << " grad.grad_inc = " << grad.grad_inc << " grad.factor " << grad.factor - // << endl; -#endif for (int y = dimY - 1; y >= 0; y--) { for (int x = 0; x < dimX; x++) { // NB: gradient lies in range 0-1 for all types, and is applied when necessary... @@ -1273,11 +1060,6 @@ void Landscape::setEnvGradient(Species* pSpecies, bool initial) break; } } -#if RSDEBUG - //DEBUGLOG << "Landscape::setEnvGradient(): y=" << y << " x=" << x - // << " dist_from_opt=" << dist_from_opt << " rasterType=" << rasterType << " hab=" << hab - // << endl; -#endif if (habK > 0.0) { // suitable cell if (initial) { // set local environmental deviation cells[y][x]->setEnvDev((float)pRandom->Random() * (2.0f) - 1.0f); @@ -1285,11 +1067,6 @@ void Landscape::setEnvGradient(Species* pSpecies, bool initial) dist_from_opt = (float)(fabs((double)grad.opt_y - (double)y)); dev = cells[y][x]->getEnvDev(); envval = 1.0 - dist_from_opt * grad.grad_inc + dev * grad.factor; -#if RSDEBUG - //DEBUGLOG << "Landscape::setEnvGradient(): y=" << y << " x=" << x - // << " dist_from_opt=" << dist_from_opt << " dev=" << dev << " p=" << p - // << endl; -#endif if (envval < 0.000001) envval = 0.0; if (envval > 1.0) envval = 1.0; } @@ -1325,11 +1102,6 @@ void Landscape::updateLocalStoch(void) { for (int x = 0; x < dimX; x++) { if (cells[y][x] != 0) { // not a no-data cell randpart = (float)(pRandom->Normal(0.0, env.std) * sqrt(1.0 - (env.ac * env.ac))); -#if RSDEBUG - //DEBUGLOG << "Landscape::updateLocalStoch(): y=" << y << " x=" << x - // << " env.std= " << env.std << " env.ac= " << env.ac << " randpart= " << randpart - // << endl; -#endif cells[y][x]->updateEps((float)env.ac, randpart); } } @@ -1364,11 +1136,6 @@ void Landscape::resetEffCosts(void) { void Landscape::setDynamicLand(bool dyn) { dynamic = dyn; } void Landscape::addLandChange(landChange c) { -#if RSDEBUG - //DebugGUI(("Landscape::addLandChange(): chgnum=" + Int2Str(c.chgnum) - // + " chgyear=" + Int2Str(c.chgyear) - // ).c_str()); -#endif landchanges.push_back(c); } @@ -1412,13 +1179,6 @@ int Landscape::readLandChange(int filenum, bool costs) simParams sim = paramsSim->getSim(); if (filenum < 0) return 19; - - //if (patchModel && (rasterType == 0 || rasterType == 2)) { - // if (filenum == 0) { // first change - // createPatchChgMatrix(); - // } - // pchseq = patchCount(); - //} if (patchModel) pchseq = patchCount(); #if !RS_RCPP @@ -1543,11 +1303,6 @@ int Landscape::readLandChange(int filenum, bool costs) } #endif } -#if RSDEBUG - //DebugGUI(("Landscape::readLandscape(): x=" + Int2Str(x) + " y=" + Int2Str(y) - // + " h=" + Int2Str(h) + " p=" + Int2Str(p) - //).c_str()); -#endif if (cells[y][x] != 0) { // not a no data cell (in initial landscape) if (h == habnodata) { // invalid no data cell in change map hfile.close(); hfile.clear(); @@ -1687,10 +1442,6 @@ int Landscape::readLandChange(int filenum, bool costs) } #endif } -#if RSDEBUG - //MemoLine(("y=" + Int2Str(y) + " x=" + Int2Str(x) + " hfloat=" + Float2Str(hfloat) - // + " p=" + Int2Str(p)).c_str()); -#endif if (cells[y][x] != 0) { // not a no data cell (in initial landscape) if (h == habnodata) { // invalid no data cell in change map hfile.close(); hfile.clear(); @@ -1795,14 +1546,6 @@ void Landscape::createPatchChgMatrix(void) } } patchChgMatrix[y][x][2] = 0; -#if RSDEBUG - //DebugGUI(("Landscape::createPatchChgMatrix(): y=" + Int2Str(y) - // + " x=" + Int2Str(x) - // + " patchChgMatrix[y][x][0]=" + Int2Str(patchChgMatrix[y][x][0]) - // + " [1]=" + Int2Str(patchChgMatrix[y][x][1]) - // + " [2]=" + Int2Str(patchChgMatrix[y][x][2]) - // ).c_str()); -#endif } } } @@ -1820,15 +1563,6 @@ void Landscape::recordPatchChanges(int landIx) { chg.oldpatch = patchChgMatrix[y][x][2]; chg.newpatch = patchChgMatrix[y][x][0]; patchchanges.push_back(chg); -#if RSDEBUG - //DebugGUI(("Landscape::recordPatchChanges(): landIx=" + Int2Str(landIx) - // + " chg.chgnum=" + Int2Str(chg.chgnum) - // + " chg.x=" + Int2Str(chg.x) - // + " chg.y=" + Int2Str(chg.y) - // + " chg.oldpatch=" + Int2Str(chg.oldpatch) - // + " chg.newpatch=" + Int2Str(chg.newpatch) - // ).c_str()); -#endif } } else { // any other change @@ -1838,15 +1572,6 @@ void Landscape::recordPatchChanges(int landIx) { chg.oldpatch = patchChgMatrix[y][x][1]; chg.newpatch = patchChgMatrix[y][x][2]; patchchanges.push_back(chg); -#if RSDEBUG - //DebugGUI(("Landscape::recordPatchChanges(): landIx=" + Int2Str(landIx) - // + " chg.chgnum=" + Int2Str(chg.chgnum) - // + " chg.x=" + Int2Str(chg.x) - // + " chg.y=" + Int2Str(chg.y) - // + " chg.oldpatch=" + Int2Str(chg.oldpatch) - // + " chg.newpatch=" + Int2Str(chg.newpatch) - // ).c_str()); -#endif } } // reset cell for next landscape change @@ -1897,14 +1622,6 @@ void Landscape::createCostsChgMatrix(void) costsChgMatrix[y][x][0] = costsChgMatrix[y][x][1] = pCell->getCost(); } costsChgMatrix[y][x][2] = 0; -#if RSDEBUG - //DebugGUI(("Landscape::createCostsChgMatrix(): y=" + Int2Str(y) - // + " x=" + Int2Str(x) - // + " costsChgMatrix[y][x][0]=" + Int2Str(costsChgMatrix[y][x][0]) - // + " [1]=" + Int2Str(costsChgMatrix[y][x][1]) - // + " [2]=" + Int2Str(costsChgMatrix[y][x][2]) - // ).c_str()); -#endif } } } @@ -1925,42 +1642,15 @@ void Landscape::recordCostChanges(int landIx) { chg.oldcost = costsChgMatrix[y][x][2]; chg.newcost = costsChgMatrix[y][x][0]; costschanges.push_back(chg); -#if RSDEBUG - //DebugGUI(("Landscape::recordCostsChanges(): landIx=" + Int2Str(landIx) - // + " chg.chgnum=" + Int2Str(chg.chgnum) - // + " chg.x=" + Int2Str(chg.x) - // + " chg.y=" + Int2Str(chg.y) - // + " chg.oldcost=" + Int2Str(chg.oldcost) - // + " chg.newcost=" + Int2Str(chg.newcost) - // ).c_str()); -#endif } } else { // any other change -#if RSDEBUG - //if (x < 20 && y == 0) { - // DEBUGLOG << "Landscape::recordCostChanges(): x=" << x << " y=" << y - // << " costsChgMatrix[y][x][0]=" << costsChgMatrix[y][x][0] - // << " costsChgMatrix[y][x][1]=" << costsChgMatrix[y][x][1] - // << " costsChgMatrix[y][x][2]=" << costsChgMatrix[y][x][2] - // << endl; - //} -#endif if (costsChgMatrix[y][x][2] != costsChgMatrix[y][x][1]) { // record change of cost for current cell chg.chgnum = landIx; chg.x = x; chg.y = y; chg.oldcost = costsChgMatrix[y][x][1]; chg.newcost = costsChgMatrix[y][x][2]; costschanges.push_back(chg); -#if RSDEBUG - //DebugGUI(("Landscape::recordCostsChanges(): landIx=" + Int2Str(landIx) - // + " chg.chgnum=" + Int2Str(chg.chgnum) - // + " chg.x=" + Int2Str(chg.x) - // + " chg.y=" + Int2Str(chg.y) - // + " chg.oldcost=" + Int2Str(chg.oldcost) - // + " chg.newcost=" + Int2Str(chg.newcost) - // ).c_str()); -#endif } } // reset cell for next landscape change @@ -2267,12 +1957,6 @@ int Landscape::readLandscape(int fileNum, string habfile, string pchfile, string return 135; } #endif - -#if RSDEBUG - //DebugGUI(("Landscape::readLandscape(): x=" + Int2Str(x) + " y=" + Int2Str(y) - // + " h=" + Int2Str(h) + " p=" + Int2Str(p) - //).c_str()); -#endif if (h == habnodata) addNewCellToLand(x, y, -1); // add cell only to landscape else { @@ -2374,10 +2058,6 @@ if (patchModel) #endif } //end if patchmodel -#if RSDEBUG -//MemoLine(("y=" + Int2Str(y) + " x=" + Int2Str(x) + " hfloat=" + Float2Str(hfloat) -// + " p=" + Int2Str(p)).c_str()); -#endif if (h == habnodata) { addNewCellToLand(x, y, -1); // add cell only to landscape } @@ -2526,10 +2206,6 @@ else { // couldn't read from hfile } #endif } -#if RSDEBUG - //MemoLine(("y=" + Int2Str(y) + " x=" + Int2Str(x) + " hfloat=" + Float2Str(hfloat) - // + " p=" + Int2Str(p)).c_str()); -#endif if (h == habnodata) { addNewCellToLand(x, y, -1); // add cell only to landscape } @@ -2631,6 +2307,7 @@ int Landscape::readCosts(string fname) wstring header; #else string header; +#endif Cell* pCell; #if !RS_RCPP simView v = paramsSim->getViews(); @@ -2644,6 +2321,7 @@ int Landscape::readCosts(string fname) #endif #endif // open cost file + // open cost file #if !RS_RCPP || RSWIN64 costs.open(fname.c_str()); #else @@ -2653,14 +2331,6 @@ int Landscape::readCosts(string fname) costs.imbue(std::locale(costs.getloc(), new std::codecvt_utf16)); } #endif - //if (!costs.is_open()) { - // MessageDlg("COSTS IS NOT OPEN!!!!!", - // mtError, TMsgDlgButtons() << mbRetry,0); - //} - //else { - // MessageDlg("Costs is open!", - // mtError, TMsgDlgButtons() << mbRetry,0); - //} // read headers and check that they correspond to the landscape ones costs >> header; #if RS_RCPP @@ -2684,6 +2354,7 @@ int Landscape::readCosts(string fname) costs >> maxXcost >> header >> maxYcost >> header >> minLongCost; costs >> header >> minLatCost >> header >> resolCost >> header >> NODATACost; + #if !RS_RCPP MemoLine("Loading costs map. Please wait..."); #endif @@ -2710,19 +2381,10 @@ int Landscape::readCosts(string fname) } #endif if (hc < 1 && hc != NODATACost) { -#if RSDEBUG -#if BATCH - // DEBUGLOG << "Landscape::readCosts(): x=" << x << " y=" << y - // << " fcost=" << fcost << " hc=" << hc - // << endl; -#endif -#endif #if RS_RCPP && !R_CMD Rcpp::Rcout << "Cost map my only contain values of 1 or higher, but found " << fcost << "." << endl; #endif // error - zero / negative cost not allowed -// MessageDlg("Error in the costs map file : zero or negative cost detected." -// , mtError, TMsgDlgButtons() << mbOK,0); costs.close(); costs.clear(); return -999; } @@ -2742,7 +2404,7 @@ if (costs.eof()) { } else EOFerrorR(fname); #else - MemoLine("Costs map loaded."); +MemoLine("Costs map loaded."); #endif costs.close(); costs.clear(); @@ -2767,40 +2429,16 @@ rasterdata CheckRasterFile(string fname) infile.open(fname.c_str()); if (infile.is_open()) { infile >> header >> r.ncols; -#if RSDEBUG - DebugGUI(("CheckRasterFile(): header=" + header + " r.ncols=" + Int2Str(r.ncols) - ).c_str()); -#endif if (header != "ncols" && header != "NCOLS") r.errors++; infile >> header >> r.nrows; -#if RSDEBUG - DebugGUI(("CheckRasterFile(): header=" + header + " r.nrows=" + Int2Str(r.nrows) - ).c_str()); -#endif if (header != "nrows" && header != "NROWS") r.errors++; infile >> header >> r.xllcorner; -#if RSDEBUG - DebugGUI(("CheckRasterFile(): header=" + header + " r.xllcorner=" + Float2Str(r.xllcorner) - ).c_str()); -#endif if (header != "xllcorner" && header != "XLLCORNER") r.errors++; infile >> header >> r.yllcorner; -#if RSDEBUG - DebugGUI(("CheckRasterFile(): header=" + header + " r.yllcorner=" + Float2Str(r.yllcorner) - ).c_str()); -#endif if (header != "yllcorner" && header != "YLLCORNER") r.errors++; infile >> header >> r.cellsize; -#if RSDEBUG - DebugGUI(("CheckRasterFile(): header=" + header + " r.cellsize=" + Int2Str(r.cellsize) - ).c_str()); -#endif if (header != "cellsize" && header != "CELLSIZE") r.errors++; infile >> header >> inint; -#if RSDEBUG - DebugGUI(("CheckRasterFile(): header=" + header + " inint=" + Int2Str(inint) - ).c_str()); -#endif if (header != "NODATA_value" && header != "NODATA_VALUE") r.errors++; infile.close(); infile.clear(); @@ -2823,9 +2461,6 @@ void Landscape::createConnectMatrix(void) { if (connectMatrix != 0) deleteConnectMatrix(); int npatches = (int)patches.size(); -#if RSDEBUG - //DEBUGLOG << "Landscape::createConnectMatrix(): npatches=" << npatches << endl; -#endif connectMatrix = new int* [npatches]; for (int i = 0; i < npatches; i++) { connectMatrix[i] = new int[npatches]; @@ -3025,10 +2660,6 @@ void Landscape::outVisits(int rep, int landNr) { outvisits << "NODATA_value -9" << endl; for (int y = dimY - 1; y >= 0; y--) { -#if RSDEBUG - //DebugGUI(("Landscape::drawLandscape(): y=" + Int2Str(y) - // + " cells[y]=" + Int2Str((int)cells[y])).c_str()); -#endif for (int x = 0; x < dimX; x++) { if (cells[y][x] == 0) { // no-data cell outvisits << "-9 "; diff --git a/src/RScore/Landscape.h b/src/RScore/Landscape.h index d73ffbd..4ef9231 100644 --- a/src/RScore/Landscape.h +++ b/src/RScore/Landscape.h @@ -140,7 +140,6 @@ class InitDist { }; - //--------------------------------------------------------------------------- struct landParams { @@ -235,6 +234,7 @@ class Landscape { void setCellArray(void); void addPatchNum(int); + std::vector getPatchNums() const { return patchnums; } void generatePatches(Species*); // create an artificial landscape void allocatePatches(Species*); // create patches for a cell-based landscape Patch* newPatch( diff --git a/src/RScore/Main.cpp b/src/RScore/Main.cpp new file mode 100644 index 0000000..42f783b --- /dev/null +++ b/src/RScore/Main.cpp @@ -0,0 +1,95 @@ +/*---------------------------------------------------------------------------- + * + * Copyright (C) 2020 Greta Bocedi, Stephen C.F. Palmer, Justin M.J. Travis, Anne-Kathleen Malchow, 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 . + * + --------------------------------------------------------------------------*/ + +#if LINUX_CLUSTER || R_CMD +#include +#else +#include +#endif + +#include +#include +#include +#include +#include "Individual.h" +#include "Community.h" +#include "RSrandom.h" +#include "Utils.h" +#include "Parameters.h" +#include "Landscape.h" +#include "Species.h" +#include "SubCommunity.h" + +using namespace std; + +void run_unit_tests() { + cout << "******* Unit test output *******" << endl; + testRSrandom(); + testIndividual(); + cout << endl << "************************" << endl; +} + +// Global vars +string habmapname, patchmapname, distnmapname; +string costmapname, genfilename; +string landFile; +paramGrad* paramsGrad; +paramStoch* paramsStoch; +paramInit* paramsInit; +paramSim* paramsSim; +RSrandom* pRandom; +ofstream DEBUGLOG; +ofstream MUTNLOG; +vector hfnames; +Species* pSpecies; +Community* pComm; +void DebugGUI(string msg) { + // nothing +} +void MemoLine(string msg) { + /// nothing +} +traitCanvas SetupTraitCanvas(void) { + traitCanvas tcanv; + for (int i = 0; i < NTRAITS; i++) { tcanv.pcanvas[i] = 0; } + return tcanv; +} +void Landscape::setLandMap(void) { } +void Landscape::drawLandscape(int rep, int yr, int landnum) { } +void Community::viewOccSuit(int year, double mn, double se) { } +void Community::draw(int rep, int yr, int gen, int landNum) { } + +#if LINUX_CLUSTER || RS_RCPP +int main(int argc, char* argv[]) +#else +int _tmain(int argc, _TCHAR* argv[]) +#endif +{ +#ifdef NDEBUG + cout << "This code is only for running tests and not meant to run in release." << endl; + return 1; +# else + assert(0.1 > 0.0); // assert does run correctly + run_unit_tests(); + cout << "All tests have completed." << endl; + return 0; // if tests succeed, we are happy +# endif // NDEBUG +} diff --git a/src/RScore/Model.cpp b/src/RScore/Model.cpp index 5a8f3e9..2d68cd5 100644 --- a/src/RScore/Model.cpp +++ b/src/RScore/Model.cpp @@ -34,11 +34,8 @@ Rcpp::List RunModel(Landscape* pLandscape, int seqsim) int RunModel(Landscape* pLandscape, int seqsim) #endif { - //int Nsuit,yr,totalInds; int yr, totalInds; - //float gradval,gradmin,gradmax; bool filesOK; - //int t0,t1; landParams ppLand = pLandscape->getLandParams(); envGradParams grad = paramsGrad->getGradient(); @@ -51,18 +48,12 @@ int RunModel(Landscape* pLandscape, int seqsim) simParams sim = paramsSim->getSim(); simView v = paramsSim->getViews(); - //t0 = time(0); - #if RSDEBUG landPix p = pLandscape->getLandPix(); DEBUGLOG << "RunModel(): reps=" << sim.reps << " ppLand.nHab=" << ppLand.nHab << " p.pix=" << p.pix << endl; - //DEBUGLOG << "RunModel(): random integers:"; - //for (int i = 0; i < 5; i++) { - // int rrrr = pRandom->IRandom(1000,2000); DEBUGLOG << " " << rrrr; - //} DEBUGLOG << endl; #endif @@ -184,17 +175,9 @@ int RunModel(Landscape* pLandscape, int seqsim) #else SubCommunity* pSubComm = pComm->addSubComm(ppp.pPatch, ppp.patchNum); // SET UP ALL SUB-COMMUNITIES #endif - // if (ppp.y >= 9995) { - // DEBUGLOG << "RunModel(): i=" << i << " pSubComm=" << pSubComm - // << endl; - // } #else pComm->addSubComm(ppp.pPatch, ppp.patchNum); // SET UP ALL SUB-COMMUNITIES #endif - // if (i == 0 || ppp.pPatch->getK() > 0.0) { - // // SET UP SUB-COMMUNITY FOR MATRIX PATCH AND ANY PATCH HAVING NON-ZERO CARRYING CAPACITY - // pComm->addSubComm(ppp.pPatch,ppp.patchNum); - // } } MemoLine("...completed..."); #if RSDEBUG @@ -320,7 +303,7 @@ int RunModel(Landscape* pLandscape, int seqsim) << " pSpecies=" << pSpecies << endl; #endif #if BATCH && RS_RCPP && !R_CMD - Rcpp::Rcout << "RunModel(): completed initialisation " << endl; + Rcpp::Rcout << "RunModel(): completed initialisation " << endl; #endif // open a new individuals file for each replicate @@ -540,7 +523,6 @@ int RunModel(Landscape* pLandscape, int seqsim) pLandscape->resetConnectMatrix(); if (ppLand.dynamic && updateland) { - // trfrRules trfr = pSpecies->getTrfr(); if (trfr.moveModel && trfr.moveType == 1) { // SMS if (!trfr.costMap) pLandscape->resetCosts(); // in case habitats have changed } @@ -601,12 +583,6 @@ int RunModel(Landscape* pLandscape, int seqsim) list_outPop.push_back(pComm->addYearToPopList(rep, yr), "rep" + std::to_string(rep) + "_year" + std::to_string(yr)); } #endif - -#if RSDEBUG - //DEBUGLOG << "RunModel(): completed RangePopOutput()" - // << " Total_Size = " << Total_Size << endl; -#endif - // apply local extinction for generation 0 only // CHANGED TO *BEFORE* RANGE & POPN OUTPUT PRODUCTION IN v1.1, // SO THAT NOS. OF JUVENILES BORN CAN BE REPORTED @@ -750,11 +726,6 @@ int RunModel(Landscape* pLandscape, int seqsim) pComm->resetPopns(); - // if (batchMode) { - // // delete the community of species using the landscape - // pComm->resetPopns(); - // } - //Reset the gradient optimum if (grad.gradient) paramsGrad->resetOptY(); @@ -878,7 +849,6 @@ int RunModel(Landscape* pLandscape, int seqsim) MemoLine("Writing final occupancy output..."); pComm->outOccupancy(); pComm->outOccSuit(v.viewGraph); - // pComm->deleteOccupancy((sim.years/sim.outInt)+1); pComm->deleteOccupancy((sim.years / sim.outIntOcc) + 1); pComm->outOccupancyHeaders(-999); MemoLine("...finished"); @@ -899,17 +869,12 @@ int RunModel(Landscape* pLandscape, int seqsim) if (sim.outInds) pComm->outInds(0, 0, 0, -999); if (sim.outputWCFstat) pComm->openWCFstatFile(pSpecies, -999); if (sim.outputPerLocusWCFstat) pComm->openWCPerLocusFstatFile(pSpecies, pLandscape, -999, 0); - if (sim.outputPairwiseFst) pComm->openPairwiseFSTFilei(pSpecies, pLandscape, -999, 0); + if (sim.outputPairwiseFst) pComm->openPairwiseFSTFile(pSpecies, pLandscape, -999, 0); MemoLine("Deleting community..."); delete pComm; pComm = 0; MemoLine("...finished"); - // Write performance data - //t1 = time(0); - //RSlog << "Simulation," << sim.simulation << "," << sim.reps << "," << sim.years - // << "," << t1-t0 << endl; - #if RS_RCPP && !R_CMD return list_outPop; #else @@ -953,7 +918,7 @@ bool CheckDirectory(void) void PreReproductionOutput(Landscape* pLand, Community* pComm, int rep, int yr, int gen) { #if RSDEBUG -landParams ppLand = pLand->getLandParams(); + landParams ppLand = pLand->getLandParams(); #endif simParams sim = paramsSim->getSim(); simView v = paramsSim->getViews(); @@ -968,10 +933,6 @@ landParams ppLand = pLand->getLandParams(); << endl; #endif -#if RSDEBUG - //DEBUGLOG << "PreReproductionOutput(): 22222 " << endl; -#endif - traitCanvas tcanv; for (int i = 0; i < NTRAITS; i++) { tcanv.pcanvas[i] = 0; @@ -989,23 +950,8 @@ landParams ppLand = pLand->getLandParams(); { pComm->outTraits(tcanv, pSpecies, rep, yr, gen); } - -#if RSDEBUG - //DEBUGLOG << "PreReproductionOutput(): 33333 " << endl; -#endif - if (sim.outOccup && yr % sim.outIntOcc == 0 && gen == 0) pComm->updateOccupancy(yr / sim.outIntOcc, rep); - -#if RSDEBUG - //DEBUGLOG << "PreReproductionOutput(): 88888 " << endl; -#endif - -// Remaining graphical output actions are performed for GUI only - -#if RSDEBUG -//DEBUGLOG << "PreReproductionOutput(): finished " << endl; -#endif } //For outputs and population visualisations pre-reproduction @@ -1164,9 +1110,7 @@ void OutParameters(Landscape* pLandscape) // << " habitat map: " << chg.habfile << endl; } } - outPar << endl << "SPECIES DISTRIBUTION LOADED: \t"; - //if (sim.initDistLoaded) if (ppLand.spDist) { outPar << "yes" << endl; @@ -1485,7 +1429,6 @@ void OutParameters(Landscape* pLandscape) else outPar << "K "; outPar << k << endl; } - emigTraits ep0, ep1; string sexdept = "SEX-DEPENDENT: "; string stgdept = "STAGE-DEPENDENT: "; @@ -1495,7 +1438,6 @@ void OutParameters(Landscape* pLandscape) outPar << endl << "DISPERSAL - EMIGRATION:\t"; if (emig.densDep) { outPar << "density-dependent" << endl; - if (emig.sexDep) { outPar << sexdept << "yes" << endl; if (emig.stgDep) { @@ -1746,7 +1688,6 @@ void OutParameters(Landscape* pLandscape) outPar << "MAX. No. OF STEPS:\t "; if (ssteps.maxSteps == 99999999) outPar << "not applied" << endl; else outPar << ssteps.maxSteps << endl; - if (sett.sexDep) { nsexes = 2; outPar << sexdept << "yes" << endl; @@ -1863,7 +1804,6 @@ void OutParameters(Landscape* pLandscape) } // Genetics - outPar << endl << "GENETICS:" << endl; set traitList = pSpecies->getTraitTypes(); diff --git a/src/RScore/NeutralStatsManager.cpp b/src/RScore/NeutralStatsManager.cpp index bcd59b5..3c0b15a 100644 --- a/src/RScore/NeutralStatsManager.cpp +++ b/src/RScore/NeutralStatsManager.cpp @@ -29,7 +29,7 @@ // ---------------------------------------------------------------------------------------- NeutralStatsManager::NeutralStatsManager(set const& patchList, const int nLoci) { - this->_fst_matrix = PatchMatrix(patchList.size(), patchList.size()); + this->_fst_matrix = PatchMatrix(static_cast(patchList.size()), static_cast(patchList.size())); globalAlleleTable.reserve(nLoci); //don't have to be pointers, not shared or moved } @@ -95,7 +95,7 @@ void NeutralStatsManager::resetGlobalAlleleTable() { void NeutralStatsManager::setLociDiversityCounter(set const& patchList, const int nInds, Species* pSpecies, Landscape* pLandscape) { - unsigned int i, j, k; + int i, j; const int nLoci = pSpecies->getNPositionsForTrait(SNP); const int nAlleles = (int)pSpecies->getSpTrait(SNP)->getMutationParameters().find(MAX)->second; const int chromosomes = (pSpecies->isDiploid() ? 2 : 1); @@ -282,17 +282,17 @@ void NeutralStatsManager::calculateFstatWC(set const& patchList, const int int patchSize = pPop->sampleSize(); if (patchSize) { extantPs++; - sum_weights += (patchSize * patchSize / nInds); + sum_weights += (patchSize * patchSize / static_cast(nInds)); } } _n_extantPopulations = extantPs; _n_individuals = nInds; - n_bar = nInds / extantPs; + n_bar = nInds / static_cast(extantPs); n_c = (nInds - sum_weights) / (extantPs - 1); - inverse_n_bar = 1 / (n_bar - 1); - inverse_n_total = 1 / nInds; + inverse_n_bar = 1.0 / (n_bar - 1); + inverse_n_total = 1.0 / nInds; double var; double s2, p_bar, h_bar; @@ -302,9 +302,9 @@ void NeutralStatsManager::calculateFstatWC(set const& patchList, const int double a = 0, b = 0, c = 0, x; - for (unsigned int thisLocus = 0; thisLocus < nLoci; ++thisLocus) { + for (int thisLocus = 0; thisLocus < nLoci; ++thisLocus) { - for (unsigned int allele = 0; allele < nAlleles; ++allele) { + for (int allele = 0; allele < nAlleles; ++allele) { s2 = p_bar = h_bar = 0; @@ -359,7 +359,7 @@ void NeutralStatsManager::calculateFstatWC_MS(set const& patchList, const i int patchSize = pPop->sampleSize(); if (patchSize) { extantPs++; - sum_weights += (patchSize * patchSize / nInds); + sum_weights += (patchSize * patchSize / static_cast(nInds)); } } @@ -371,20 +371,20 @@ void NeutralStatsManager::calculateFstatWC_MS(set const& patchList, const i //p = _alleleFreqTable //pb = _globalAlleleFreq - vector alploc(nLoci); + vector alploc(nLoci); unsigned int** alploc_table = new unsigned int* [nLoci]; - for (unsigned int i = 0; i < nLoci; ++i) + for (int i = 0; i < nLoci; ++i) alploc_table[i] = new unsigned int[nAlleles]; - unsigned int tot_num_allele = 0; + int tot_num_allele = 0; - for (unsigned int l = 0; l < nLoci; ++l) { + for (int l = 0; l < nLoci; ++l) { alploc[l] = 0; - for (unsigned int cnt, a = 0; a < nAlleles; ++a) { + for (int cnt, a = 0; a < nAlleles; ++a) { cnt = 0; @@ -411,13 +411,13 @@ void NeutralStatsManager::calculateFstatWC_MS(set const& patchList, const i vector SSP(tot_num_allele); vector SSi(tot_num_allele); - unsigned int all_cntr = 0; + int all_cntr = 0; double het, freq, var; - for (unsigned int l = 0; l < nLoci; ++l) { + for (int l = 0; l < nLoci; ++l) { - for (unsigned int a = 0; a < nAlleles && all_cntr < tot_num_allele; ++a) { + for (int a = 0; a < nAlleles && all_cntr < tot_num_allele; ++a) { if (alploc_table[l][a] == 0) continue; //do not consider alleles not present in the pop @@ -475,7 +475,7 @@ void NeutralStatsManager::calculateFstatWC_MS(set const& patchList, const i _fit_WC_loc.resize(nLoci); if (tot_num_allele != nLoci) { - for (unsigned int i = 0; i < tot_num_allele; ++i) { + for (int i = 0; i < tot_num_allele; ++i) { MSG[i] = SSG[i] / (2 * nInds); sigw[i] = MSG[i]; //wasted! @@ -499,11 +499,11 @@ void NeutralStatsManager::calculateFstatWC_MS(set const& patchList, const i // cout<<" computing sigma per locus\n"; - for (unsigned int allcntr = 0, i = 0; i < nLoci; ++i) { + for (int allcntr = 0, i = 0; i < nLoci; ++i) { lsiga = lsigb = lsigw = 0; - for (unsigned int l = 0; l < alploc[i]; ++l) { + for (int l = 0; l < alploc[i]; ++l) { lsiga += siga[allcntr]; lsigb += sigb[allcntr]; lsigw += sigw[allcntr]; @@ -529,7 +529,7 @@ void NeutralStatsManager::calculateFstatWC_MS(set const& patchList, const i _fis_WC = 0; } - for (unsigned int i = 0; i < nLoci; ++i) + for (int i = 0; i < nLoci; ++i) delete[]alploc_table[i]; delete[]alploc_table; } @@ -550,7 +550,7 @@ void NeutralStatsManager::setFstMatrix(set const& patchList, const int nInd copy(patchList.begin(), patchList.end(), std::back_inserter(patchVect)); //needs to be in vector to iterate over, copy preserves order - int nPatches = patchList.size(); + int nPatches = static_cast(patchList.size()); //initialise @@ -564,7 +564,7 @@ void NeutralStatsManager::setFstMatrix(set const& patchList, const int nInd double* pop_weights = new double[nPatches]; double* pop_sizes = new double[nPatches]; double** numerator = new double* [nPatches]; - for (unsigned int i = 0; i < nPatches; i++) numerator[i] = new double[nPatches]; + for (int i = 0; i < nPatches; i++) numerator[i] = new double[nPatches]; double tot_size; double numerator_W = 0; double denominator = 0; @@ -572,20 +572,20 @@ void NeutralStatsManager::setFstMatrix(set const& patchList, const int nInd tot_size = nInds * 2; //diploid - for (unsigned int i = 0; i < nPatches; ++i) { + for (int i = 0; i < nPatches; ++i) { const auto patch = pLandscape->findPatch(patchVect[i]); const auto pPop = (Population*)patch->getPopn((intptr)pSpecies); pop_sizes[i] = pPop->sampleSize() * 2; pop_weights[i] = pop_sizes[i] - (pop_sizes[i] * pop_sizes[i] / tot_size); //n_ic in Weir & Hill 2002 sum_weights += pop_weights[i]; - for (unsigned int j = 0; j < nPatches; j++) + for (int j = 0; j < nPatches; j++) numerator[i][j] = 0; } double p, pq, var, num; - for (unsigned int i = 0; i < nPatches; ++i) { + for (int i = 0; i < nPatches; ++i) { if (!pop_sizes[i]) continue; @@ -598,9 +598,9 @@ void NeutralStatsManager::setFstMatrix(set const& patchList, const int nInd // // << endl; //#endif - for (unsigned int l = 0; l < nLoci; ++l) { + for (int l = 0; l < nLoci; ++l) { - for (unsigned int u = 0; u < nAlleles; ++u) { + for (int u = 0; u < nAlleles; ++u) { p = pPop->getAlleleFrequency(l, u); //p_liu @@ -630,7 +630,7 @@ void NeutralStatsManager::setFstMatrix(set const& patchList, const int nInd }// end for locus }//end for pop - for (unsigned int i = 0; i < nPatches; ++i) { + for (int i = 0; i < nPatches; ++i) { if (!pop_sizes[i]) continue; if(denominator != 0) _fst_matrix.set(i, i, 1 - (numerator[i][i] * sum_weights / denominator)); @@ -646,15 +646,15 @@ void NeutralStatsManager::setFstMatrix(set const& patchList, const int nInd //pairwise Fst: double pi, pj; - for (unsigned int l = 0; l < nLoci; ++l) - for (unsigned int u = 0; u < nAlleles; ++u) - for (unsigned int i = 0; i < nPatches - 1; ++i) { + for (int l = 0; l < nLoci; ++l) + for (int u = 0; u < nAlleles; ++u) + for (int i = 0; i < nPatches - 1; ++i) { if (!pop_sizes[i]) continue; const auto patch = pLandscape->findPatch(patchVect[i]); const auto pPopI = (Population*)patch->getPopn((intptr)pSpecies); - for (unsigned int j = i + 1; j < nPatches; ++j) { + for (int j = i + 1; j < nPatches; ++j) { if (!pop_sizes[j]) continue; const auto patch = pLandscape->findPatch(patchVect[j]); const auto pPopJ = (Population*)patch->getPopn((intptr)pSpecies); @@ -665,9 +665,9 @@ void NeutralStatsManager::setFstMatrix(set const& patchList, const int nInd } } - for (unsigned int i = 0; i < nPatches - 1; ++i) { + for (int i = 0; i < nPatches - 1; ++i) { if (!pop_sizes[i]) continue; - for (unsigned int j = i + 1; j < nPatches; ++j) { + for (int j = i + 1; j < nPatches; ++j) { if (!pop_sizes[j]) continue; if (denominator != 0) _fst_matrix.set(i, j, 1 - ((numerator[i][j] * sum_weights) / (2 * denominator))); @@ -684,7 +684,7 @@ void NeutralStatsManager::setFstMatrix(set const& patchList, const int nInd } delete[] pop_weights; delete[] pop_sizes; - for (unsigned int i = 0; i < nPatches; i++) delete[] numerator[i]; + for (int i = 0; i < nPatches; i++) delete[] numerator[i]; delete[] numerator; } diff --git a/src/RScore/NeutralStatsManager.h b/src/RScore/NeutralStatsManager.h index 23ce1bb..4b269e0 100644 --- a/src/RScore/NeutralStatsManager.h +++ b/src/RScore/NeutralStatsManager.h @@ -15,9 +15,7 @@ struct PatchMatrix public: - PatchMatrix() {}; - - PatchMatrix(unsigned int rows, unsigned int cols) : _rows(0), _cols(0), _length(0), _val(0) { + PatchMatrix(int rows = 0, int cols = 0) : _rows(0), _cols(0), _length(0), _val(0) { _length = rows * cols; _val.resize(_length); _rows = rows; _cols = cols; diff --git a/src/RScore/Parameters.cpp b/src/RScore/Parameters.cpp index 37d9b5f..6e25de6 100644 --- a/src/RScore/Parameters.cpp +++ b/src/RScore/Parameters.cpp @@ -186,12 +186,6 @@ float paramInit::getProp(short stg) { } void paramInit::addInitInd(initInd iind) { -#if RSDEBUG - //DebugGUI(("paramInit::addInitInd(): iind.patchID=" + Int2Str(iind.patchID) - // + " iind.x=" + Int2Str(iind.x) - // + " iind.y=" + Int2Str(iind.y) - // ).c_str()); -#endif initinds.push_back(iind); } @@ -401,6 +395,104 @@ string paramSim::getDir(int option) { return s; } +const sex_t stringToSex(const std::string& str) { + if (str == "female") return FEM; + else if (str == "male") return MAL; + else throw logic_error("Traits file: ERROR - sex can either be 'female' or 'male'."); +} + +set convertStringToPatches(const string& str, const int& nb_rnd_patches, const vector& existingPatches) { + + set patches; + + if (str == "random") { + if (nb_rnd_patches > existingPatches.size()) { + throw logic_error("Genetics file: ERROR - Nb of patches to randomly sample exceeds nb of existing patches."); + } else { + // Sample without replacement + std::sample( + existingPatches.begin(), + existingPatches.end(), + std::inserter(patches, patches.end()), + nb_rnd_patches, + pRandom->getRNG() + ); + } + } else if (str == "all") { + // Copy all patches into sampled patches + std::copy(existingPatches.begin(), + existingPatches.end(), + std::inserter(patches, patches.end()) + ); + } else { + // comma-separated list of patches + stringstream ss(str); + string strPch; + int pch; + bool patchExists; + // Read comma-separated values + while (std::getline(ss, strPch, ',')) { + pch = std::stoi(strPch); + patchExists = std::find(existingPatches.begin(), existingPatches.end(), pch) != existingPatches.end(); + if (!patchExists) + throw logic_error("Genetics file: ERROR - sampled patch does not exist."); + else { + patches.insert(pch); + } + } + } + return patches; +} + +set convertStringToChromosomeEnds(string str, int genomeSize) { + set chromosomeEnds; + if (str == "#") + chromosomeEnds.insert(genomeSize - 1); // last position in genome + else { + // Parse comma-separated list from input string + stringstream ss(str); + string strPos; + int pos; + // Read comma-separated positions + while (std::getline(ss, strPos, ',')) { + pos = std::stoi(strPos); + if (pos > genomeSize) + throw logic_error("Genetics file: ERROR - chromosome ends must not exceed genome size."); + else { + chromosomeEnds.insert(pos); + } + } + } + return chromosomeEnds; +} + +set convertStringToStages(const string& str, const int& nbStages) { + set stages; + if (str == "all") { + for (int stg = 0; stg < nbStages; ++stg) { + stages.insert(stg); + } + } + else { + // Parse comma-separated list from input string + stringstream ss(str); + string strStg; + int stg; + // Read comma-separated values + while (std::getline(ss, strStg, ',')) { + stg = std::stoi(strStg); + if (stg > nbStages - 1) + throw logic_error("Genetics file: ERROR - sampled stage exceeds number of stages."); + else { + stages.insert(stg); + } + } + } + return stages; +} + + + #if RS_RCPP bool paramSim::getReturnPopRaster(void) { return ReturnPopRaster; } bool paramSim::getCreatePopFile(void) { return CreatePopFile; } diff --git a/src/RScore/Parameters.h b/src/RScore/Parameters.h index 23d8bda..344f8ac 100644 --- a/src/RScore/Parameters.h +++ b/src/RScore/Parameters.h @@ -64,8 +64,6 @@ using namespace std; #include "RSrandom.h" -class Landscape; - #define NODATACOST 100000 // cost to use in place of nodata value for SMS; #define ABSNODATACOST 100 // cost to use in place of nodata value for SMS; // when boundaries are absorbing @@ -159,11 +157,8 @@ typedef enum { KERNEL, SMS, CRW} movement_t; //GeneType convertToGeneType(const string& ); -float convertParameters(string, string); -bool iequals(std::string_view lhs, std::string_view rhs); -set convertStringToSet(string); -set convertStringToPatches(string, int, Landscape*); -set convertStringToStages(string); +set convertStringToPatches(const string&, const int&, const vector&); +set convertStringToStages(const string&, const int&); set convertStringToChromosomeEnds(string, int); //sex types @@ -451,5 +446,7 @@ extern ofstream DEBUGLOG; void DebugGUI(string); #endif +extern RSrandom* pRandom; + //--------------------------------------------------------------------------- #endif diff --git a/src/RScore/Population.cpp b/src/RScore/Population.cpp index 751c465..aec1141 100644 --- a/src/RScore/Population.cpp +++ b/src/RScore/Population.cpp @@ -40,10 +40,6 @@ Population::Population(void) { Population::Population(Species* pSp, Patch* pPch, int ninds, int resol) { // constructor for a Population of a specified size -#if RSDEBUG -//DEBUGLOG << "Population::Population(): this=" << this -// << " pPch=" << pPch << " ninds="<< ninds << endl; -#endif int n, nindivs, age = 0, minage, maxage, nAges = 0; int cumtotal = 0; @@ -60,7 +56,7 @@ Population::Population(Species* pSp, Patch* pPch, int ninds, int resol) pSpecies = pSp; pPatch = pPch; // record the new population in the patch - patchPopn pp; + patchPopn pp = patchPopn(); pp.pSp = (intptr)pSpecies; pp.pPop = (intptr)this; pPatch->addPopn(pp); #if RSDEBUG @@ -233,7 +229,7 @@ Population::~Population(void) { traitsums Population::getTraits(Species* pSpecies) { int g; - traitsums ts; + traitsums ts = traitsums(); for (int i = 0; i < NSEXES; i++) { ts.ninds[i] = 0; ts.sumD0[i] = ts.ssqD0[i] = 0.0; @@ -382,7 +378,7 @@ void Population::updateAlleleTable() { std::for_each(alleleTable.begin(), alleleTable.end(), [&](NeutralData& v) -> void { - v.setFrequencies(sampledInds.size() * 2); + v.setFrequencies(static_cast(sampledInds.size()) * 2); //v->divideHeteros(sampledInds.size()); //weir and cockerham doesn't need this division?? }); } @@ -460,10 +456,10 @@ double Population::computeHs() { popStats Population::getStats(void) { - popStats p; + popStats p = popStats(); int ninds; float fec; - bool breeders[2]; breeders[0] = breeders[1] = false; + bool breeders[2] = { false, false }; demogrParams dem = pSpecies->getDemogr(); p.pSpecies = pSpecies; p.pPatch = pPatch; @@ -522,12 +518,6 @@ int Population::totalPop(void) { for (int stg = 0; stg < nStages; stg++) { for (int sex = 0; sex < nSexes; sex++) { t += nInds[stg][sex]; -#if RSDEBUG - //DEBUGLOG << "Population::totalPop(): this=" << this - // << " stg=" << stg << " sex=" << sex - // << " nInds[stg][sex]=" << nInds[stg][sex] << " t=" << t - // << endl; -#endif } } return t; @@ -729,11 +719,6 @@ void Population::reproduction(const float localK, const float envval, const int case 0: // asexual model for (int i = 0; i < ninds; i++) { stage = inds[i]->breedingFem(); -#if RSDEBUG - //DEBUGLOG << "Population::reproduction():" - // << " i=" << i << " ID=" << inds[i]->getId() << " stage=" << stage - // << endl; -#endif if (stage > 0) { // female of breeding age if (dem.stageStruct) { // determine whether she must miss current breeding attempt @@ -869,10 +854,6 @@ void Population::reproduction(const float localK, const float envval, const int int rrr = 0; if (nmales > 1) rrr = pRandom->IRandom(0, nmales - 1); father = fathers[rrr]; -#if RSDEBUG - //DEBUGLOG << "Population::reproduction(): i = " << i << " rrr = " << rrr - // << " father = " << father << endl; -#endif pCell = pPatch->getRandomCell(); for (int j = 0; j < njuvs; j++) { @@ -910,20 +891,8 @@ void Population::reproduction(const float localK, const float envval, const int } // end of switch (dem.repType) -#if RSDEBUG -//DEBUGLOG << "Population::reproduction(): before reprodn. " << " inds.size() = " << inds.size() -// << endl; -#endif - // THIS MAY NOT BE CORRECT FOR MULTIPLE SPECIES IF THERE IS SOME FORM OF // CROSS-SPECIES DENSITY-DEPENDENT FECUNDITY - -#if RSDEBUG -//DEBUGLOG << "Population::reproduction(): after reprodn. this = " << this -// << " juvs.size() = " << juvs.size() << " inds.size() = " << inds.size() -// << endl; -#endif - } // Following reproduction of ALL species, add juveniles to the population prior to dispersal @@ -957,7 +926,7 @@ void Population::fledge(void) } Individual* Population::sampleInd() const { - int index = pRandom->IRandom(0, inds.size() - 1); + int index = pRandom->IRandom(0, static_cast(inds.size() - 1)); return inds[index]; } @@ -984,7 +953,7 @@ void Population::sampleIndsWithoutReplacement(string n, const set& sampleSt } int Population::sampleSize() const { - return sampledInds.size(); + return static_cast(sampledInds.size()); } set Population::getIndividualsInStage(int stage) { @@ -1160,10 +1129,7 @@ void Population::emigration(float localK) } // end of no individual variability disp = pRandom->Bernoulli(Pdisp); -#if RSDEBUG - //DEBUGLOG << "Population::emigration(): i=" << i << " sex=" << ind.sex << " stage=" << ind.stage - // << " Pdisp=" << Pdisp << " disp=" << disp << endl; -#endif + if (disp == 1) { // emigrant inds[i]->setStatus(1); } @@ -1181,15 +1147,8 @@ void Population::allEmigrate(void) { // If an Individual has been identified as an emigrant, remove it from the Population disperser Population::extractDisperser(int ix) { - disperser d; + disperser d = disperser(); indStats ind = inds[ix]->getStats(); -#if RSDEBUG - //if (ind.status > 0) { - // DEBUGLOG << "Population::extractDisperser(): ix = " << ix << " inds[ix] = " << inds[ix] - // << " indId = " << inds[ix]->getId() << " ind.status = " << ind.status - // << endl; - //} -#endif if (ind.status == 1) { // emigrant d.pInd = inds[ix]; d.yes = true; inds[ix] = 0; @@ -1231,17 +1190,12 @@ void Population::addSettleTraitsForInd(int ix, settleTraits& avgSettleTraits) { // if it is a settler, return its new location and remove it from the current population // otherwise, leave it in the matrix population for possible reporting before deletion disperser Population::extractSettler(int ix) { - disperser d; + disperser d = disperser(); Cell* pCell; - //Patch* pPatch; +//Patch* pPatch; indStats ind = inds[ix]->getStats(); - pCell = inds[ix]->getLocn(1); -#if RSDEBUG - //DEBUGLOG << "Population::extractSettler(): ix = " << ix << " inds[ix] = " << inds[ix] - // << " indId = " << inds[ix]->getId() << " ind.status = " << ind.status << endl; -#endif d.pInd = inds[ix]; d.pCell = pCell; d.yes = false; if (ind.status == 4 || ind.status == 5) { // settled d.yes = true; @@ -1252,17 +1206,10 @@ disperser Population::extractSettler(int ix) { } // Add a specified individual to the new/current dispersal group -// Add a specified individual to the population void Population::recruit(Individual* pInd) { inds.push_back(pInd); indStats ind = pInd->getStats(); nInds[ind.stage][ind.sex]++; -#if RSDEBUG - //DEBUGLOG << "Population::recruit(): patchNum=" << pPatch->getPatchNum() - // << " indId=" << pInd->getId() - // << " nInds[" << ind.stage << "][" << ind.sex << "]=" << nInds[ind.stage][ind.sex] - // << endl; -#endif } //--------------------------------------------------------------------------- @@ -1285,7 +1232,8 @@ int Population::transfer(Landscape* pLandscape, short landIx) Cell* pCell = 0; indStats ind; Population* pNewPopn = 0; - locn newloc, nbrloc; + locn newloc = locn(); + locn nbrloc = locn(); landData ppLand = pLandscape->getLandData(); short reptype = pSpecies->getRepType(); @@ -1295,7 +1243,6 @@ int Population::transfer(Landscape* pLandscape, short landIx) settleTraits settDD; settlePatch settle; simParams sim = paramsSim->getSim(); - // each individual takes one step // for dispersal by kernel, this should be the only step taken int ninds = (int)inds.size(); @@ -1349,7 +1296,6 @@ int Population::transfer(Landscape* pLandscape, short landIx) // << " ndispersers=" << ndispersers << endl; #endif -// each individual which has reached a potential patch decides whether to settle for (int i = 0; i < ninds; i++) { ind = inds[i]->getStats(); if (ind.sex == 0) othersex = 1; else othersex = 0; @@ -1399,7 +1345,6 @@ int Population::transfer(Landscape* pLandscape, short landIx) else { // no requirement to find a mate mateOK = true; } - densdepOK = false; settle = inds[i]->getSettPatch(); if (sett.densDep) @@ -1508,7 +1453,6 @@ int Population::transfer(Landscape* pLandscape, short landIx) } } #endif - if (!trfr.moveModel && sett.go2nbrLocn && (ind.status == 3 || ind.status == 6)) { // for kernel-based transfer only ... @@ -1622,17 +1566,6 @@ void Population::survival0(float localK, short option0, short option1) // option1: 0 - development only (when survival is annual) // 1 - development and survival // 2 - survival only (when survival is annual) - -#if RSDEBUG -//DEBUGLOG << "Population::survival0():" -// << " pSpecies=" << pSpecies << " this=" << this << " PatchNum=" << pPatch->getPatchNum() -//#if SEASONAL -// << " season=" << season -//#endif // SEASONAL -// << " localK=" << localK << " option=" << option -// << endl; -#endif - densDepParams ddparams = pSpecies->getDensDep(); demogrParams dem = pSpecies->getDemogr(); stageParams sstruct = pSpecies->getStage(); @@ -1640,7 +1573,6 @@ void Population::survival0(float localK, short option0, short option1) // get surrent population size int ninds = (int)inds.size(); if (ninds == 0) return; - // set up local copies of species development and survival tables int nsexes; if (dem.repType == 0) nsexes = 1; else nsexes = 2; @@ -1680,22 +1612,11 @@ void Population::survival0(float localK, short option0, short option1) #endif } } - if (dem.stageStruct) { -#if RSDEBUG - // DEBUGLOG << "Population::survival0(): 2222 " - // << " ninds=" << ninds << " localK=" << localK - // << " effect of density dependence:" << endl; -#endif - // apply density dependence in development and/or survival probabilities for (int stg = 0; stg < nStages; stg++) { for (int sex = 0; sex < nsexes; sex++) { if (option1 != 2 && sstruct.devDens && stg > 0) { -#if RSDEBUG - // DEBUGLOG << "DD in DEVELOPMENT for stg=" << stg << " sex=" << sex << endl; -#endif // NB DD in development does NOT apply to juveniles, - // which must develop to stage 1 if they survive float effect = 0.0; if (sstruct.devStageDens) { // stage-specific density dependence // NOTE: matrix entries represent effect of ROW on COLUMN @@ -1777,7 +1698,6 @@ void Population::survival0(float localK, short option0, short option1) } } } - // identify which individuals die or develop #if RSDEBUG //DEBUGLOG << "Population::survival0():" << " ninds " << ninds @@ -1797,6 +1717,9 @@ void Population::survival0(float localK, short option0, short option1) if ((ind.stage == 0 && option0 < 2) || (ind.stage > 0 && option0 > 0)) { // condition for processing the stage is met... if (ind.status < 6) { // not already doomed + if (ind.sex < sex_t::FEM || ind.sex > sex_t::MAL) + // ?? MSVC believes it's important to bound check ind.sex + throw runtime_error("Individual sex is out of bounds"); double probsurv = surv[ind.stage][ind.sex]; // does the individual survive? if (pRandom->Bernoulli(probsurv)) { // survives @@ -1845,7 +1768,6 @@ void Population::survival0(float localK, short option0, short option1) // Apply survival changes to the population void Population::survival1(void) { - int ninds = (int)inds.size(); #if RSDEBUG //DEBUGLOG << "Population::survival1(): this=" << this @@ -1880,29 +1802,14 @@ void Population::survival1(void) // << endl; #endif -#if RSDEBUG -//for (int i = 0; i < inds.size(); i++) { -//DEBUGLOG << "Population::survival1():" << " inds[" << i << "] = " << inds[i] << endl; -//} -#endif -// remove pointers to dead individuals clean(); -#if RSDEBUG - //DEBUGLOG << "Population::survival1(): this=" << this - // << " patchNum=" << pPatch->getPatchNum() << " finished" - // << endl; -#endif } void Population::ageIncrement(void) { int ninds = (int)inds.size(); stageParams sstruct = pSpecies->getStage(); -#if RSDEBUG - //DEBUGLOG << "Population::ageIncrement():" << " inds.size() = " << inds.size() - // << endl; -#endif for (int i = 0; i < ninds; i++) { inds[i]->ageIncrement(sstruct.maxAge); } @@ -1914,12 +1821,6 @@ void Population::clean(void) { int ninds = (int)inds.size(); if (ninds > 0) { - // sort (inds.begin(), inds.end()); - // reverse (inds.begin(), inds.end()); - // - // while (inds.size() > 0 && inds[inds.size()-1] == NULL ) { - // inds.pop_back(); - // } // ALTERNATIVE METHOD: AVOIDS SLOW SORTING OF POPULATION std::vector survivors; // all surviving individuals for (int i = 0; i < ninds; i++) { @@ -1952,10 +1853,7 @@ bool Population::outPopHeaders(int landNr, bool patchModel) { outPop.clear(); return true; } - string name; - //landParams ppLand = pLandscape->getLandParams(); - //envStochParams env = paramsStoch->getStoch(); simParams sim = paramsSim->getSim(); envGradParams grad = paramsGrad->getGradient(); @@ -1963,7 +1861,6 @@ bool Population::outPopHeaders(int landNr, bool patchModel) { // ATTRIBUTES OF *ALL* SPECIES AS DETECTED AT MODEL LEVEL demogrParams dem = pSpecies->getDemogr(); stageParams sstruct = pSpecies->getStage(); - if (sim.batchMode) { name = paramsSim->getDir(2) + "Batch" + Int2Str(sim.batchNum) + "_" @@ -1982,13 +1879,6 @@ bool Population::outPopHeaders(int landNr, bool patchModel) { if (paramsStoch->envStoch()) writeEnv = true; if (writeEnv) outPop << "\tEpsilon\tGradient\tLocal_K"; outPop << "\tSpecies\tNInd"; -#if RSDEBUG - //DEBUGLOG << "Population::outPopHeaders(): this=" << this - // << " patchNum=" << pPatch->getPatchNum() - // << " totalPop()=" << totalPop() - // << " nStages=" << nStages << " nSexes=" << nSexes - // << endl; -#endif if (dem.stageStruct) { if (dem.repType == 0) { @@ -2005,7 +1895,6 @@ bool Population::outPopHeaders(int landNr, bool patchModel) { if (dem.repType != 0) outPop << "\tNfemales\tNmales"; } outPop << endl; - return outPop.is_open(); } @@ -2016,18 +1905,10 @@ void Population::outPopulation(int rep, int yr, int gen, float eps, { Cell* pCell; -#if RSDEBUG - //DEBUGLOG << "Population::outPopulations(): this=" << this - // << " writeEnv " << (int)writeEnv - // << endl; -#endif - // 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(); popStats p; - outPop << rep << "\t" << yr << "\t" << gen; if (patchModel) { outPop << "\t" << pPatch->getPatchNum(); @@ -2172,7 +2053,6 @@ void Population::outIndsHeaders(int rep, int landNr, bool patchModel) + "_Rep" + Int2Str(rep) + "_Inds.txt"; } outInds.open(name.c_str()); - outInds << "Rep\tYear\tRepSeason\tSpecies\tIndID\tStatus"; if (patchModel) outInds << "\tNatal_patch\tPatchID"; else outInds << "\tNatal_X\tNatal_Y\tX\tY"; @@ -2219,9 +2099,7 @@ void Population::outIndividual(Landscape* pLandscape, int rep, int yr, int gen, bool writeInd; pathSteps steps; Cell* pCell; - landParams ppLand = pLandscape->getLandParams(); - //landOrigin lim = pLandscape->getOrigin(); demogrParams dem = pSpecies->getDemogr(); emigRules emig = pSpecies->getEmig(); trfrRules trfr = pSpecies->getTrfr(); @@ -2229,7 +2107,6 @@ void Population::outIndividual(Landscape* pLandscape, int rep, int yr, int gen, short spNum = pSpecies->getSpNum(); int ninds = (int)inds.size(); - for (int i = 0; i < ninds; i++) { indStats ind = inds[i]->getStats(); if (yr == -1) { // write all initialised individuals @@ -2256,7 +2133,7 @@ void Population::outIndividual(Landscape* pLandscape, int rep, int yr, int gen, outInds << "\t" << ind.status; } pCell = inds[i]->getLocn(1); - locn loc; + locn loc = locn(); if (pCell == 0) loc.x = loc.y = -1; // beyond boundary or in no-data cell else loc = pCell->getLocn(); pCell = inds[i]->getLocn(0); @@ -2284,7 +2161,7 @@ void Population::outIndividual(Landscape* pLandscape, int rep, int yr, int gen, if (dem.stageStruct) outInds << "\t" << ind.age << "\t" << ind.stage; if (pSpecies->getNumberOfAdaptiveTraits() > 0) outInds << "\t" << inds[i]->getFitness(); - + if (emig.indVar) { emigTraits e = inds[i]->getEmigTraits(); if (emig.densDep) { @@ -2294,7 +2171,6 @@ void Population::outIndividual(Landscape* pLandscape, int rep, int yr, int gen, outInds << "\t" << e.d0; } } // end of if (emig.indVar) - if (trfr.indVar) { if (trfr.moveModel) { if (trfr.moveType == 1) { // SMS @@ -2354,7 +2230,6 @@ void Population::outIndividual(Landscape* pLandscape, int rep, int yr, int gen, } - //--------------------------------------------------------------------------- //--------------------------------------------------------------------------- //--------------------------------------------------------------------------- diff --git a/src/RScore/Population.h b/src/RScore/Population.h index c7743b1..dae449a 100644 --- a/src/RScore/Population.h +++ b/src/RScore/Population.h @@ -50,9 +50,6 @@ #include #include - //#include - //#include - //#include using namespace std; #include "Parameters.h" @@ -238,7 +235,6 @@ class Population { int countHeterozygoteLoci(); vector countLociHeterozyotes(); double computeHs(); - void updateHeteroTable(); private: short nStages; diff --git a/src/RScore/QTLTrait.cpp b/src/RScore/QTLTrait.cpp index 6c9d451..1b38fee 100644 --- a/src/RScore/QTLTrait.cpp +++ b/src/RScore/QTLTrait.cpp @@ -20,10 +20,10 @@ QTLTrait::QTLTrait(SpeciesTrait* P) switch (mutationDistribution) { case UNIFORM: { - if (!mutationParameters.count(MAX)) + if (mutationParameters.count(MAX) != 1) cout << endl << ("Error:: mutation uniform qtl distribution parameter must contain max value (e.g. max= ) \n"); - if (!mutationParameters.count(MIN)) + if (mutationParameters.count(MIN) != 1) cout << endl << ("Error:: mutation uniform qtl distribution parameter must contain min value (e.g. min= ) \n"); _mutate_func_ptr = &QTLTrait::mutateUniform; @@ -31,10 +31,10 @@ QTLTrait::QTLTrait(SpeciesTrait* P) } case NORMAL: { - if (!mutationParameters.count(MEAN)) + if (mutationParameters.count(MEAN) != 1) cout << endl << ("Error:: qtl mutation distribution set to normal so parameters must contain mean value (e.g. mean= ) \n"); - if (!mutationParameters.count(SDEV)) + if (mutationParameters.count(SDEV) != 1) cout << endl << ("Error::qtl mutation distribution set to normal so parameters must contain sdev value (e.g. sdev= ) \n"); _mutate_func_ptr = &QTLTrait::mutateNormal; @@ -55,10 +55,10 @@ QTLTrait::QTLTrait(SpeciesTrait* P) switch (initialDistribution) { case UNIFORM: { - if (!initialParameters.count(MAX)) + if (initialParameters.count(MAX) != 1) cout << endl << ("Error:: initial uniform qtl distribution parameter must contain max value (e.g. max= ) \n"); - if (!initialParameters.count(MIN)) + if (initialParameters.count(MIN) != 1) cout << endl << ("Error:: initial uniform qtl distribution parameter must contain min value (e.g. min= ) \n"); float maxD = initialParameters.find(MAX)->second; @@ -69,10 +69,10 @@ QTLTrait::QTLTrait(SpeciesTrait* P) } case NORMAL: { - if (!initialParameters.count(MEAN)) + if (initialParameters.count(MEAN) != 1) cout << endl << ("Error:: initial normal qtl distribution parameter must contain mean value (e.g. mean= ) \n"); - if (!initialParameters.count(SDEV)) + if (initialParameters.count(SDEV) != 1) cout << endl << ("Error:: initial normal qtl distribution parameter must contain sdev value (e.g. sdev= ) \n"); float mean = initialParameters.find(MEAN)->second; @@ -209,18 +209,18 @@ void QTLTrait::inheritDiploid(sex_t whichChromosome, mapfirst); - unsigned int nextBreakpoint = *it; + int nextBreakpoint = *it; auto distance = std::distance(recomPositions.begin(), it); if (distance % 2 != 0) - parentChromosome = !parentChromosome; //switch chromosome + parentChromosome = 1 - parentChromosome; //switch chromosome for (auto const& [locus, allelePair] : parentGenes) { while (locus > nextBreakpoint) { std::advance(it, 1); nextBreakpoint = *it; - parentChromosome = !parentChromosome; //switch chromosome + parentChromosome = 1 - parentChromosome; //switch chromosome } if (locus <= nextBreakpoint) { @@ -258,10 +258,10 @@ void QTLTrait::inheritInitialParameters(sex_t whichChromosome, mapsecond; @@ -273,10 +273,10 @@ void QTLTrait::inheritInitialParameters(sex_t whichChromosome, mapsecond; @@ -303,7 +303,7 @@ void QTLTrait::inheritInitialParameters(sex_t whichChromosome, mapgetPositions(); + const set positions = pSpeciesTrait->getPositions(); short ploidy = pSpeciesTrait->getPloidy(); for (auto position : positions) { @@ -318,7 +318,7 @@ void QTLTrait::initialiseNormal(float mean, float sd) { void QTLTrait::initialiseUniform(float min, float max) { - const auto positions = pSpeciesTrait->getPositions(); + const set positions = pSpeciesTrait->getPositions(); short ploidy = pSpeciesTrait->getPloidy(); for (auto position : positions) { @@ -341,7 +341,7 @@ float QTLTrait::expressAdditive() { for (auto const& [locus, allelePair] : genes) { - for (auto m : allelePair) + for (const std::shared_ptr m : allelePair) phenotype += m->getAlleleValue(); } return phenotype; diff --git a/src/RScore/QTLTrait.h b/src/RScore/QTLTrait.h index 8840abe..2a62de8 100644 --- a/src/RScore/QTLTrait.h +++ b/src/RScore/QTLTrait.h @@ -13,7 +13,7 @@ class QTLTrait : public TTrait { private: - constexpr double QTLDominanceFactor = 1.0; // that is, no dominance + const double QTLDominanceFactor = 1.0; // that is, no dominance SpeciesTrait* pSpeciesTrait; // would be better as const so immutable, but means passing positions list is heavy and can't be passed by reference diff --git a/src/RScore/README.md b/src/RScore/README.md index 10a8386..f8dcbb9 100644 --- a/src/RScore/README.md +++ b/src/RScore/README.md @@ -1,3 +1,72 @@ -# Rangeshifter core code +# RangeShifter core code -This repo contains the core simulation code for RangeShifter v2.0 +This repo contains the core simulation code for RangeShifter v2.0 and is not meant to be compiled or run on its own. + + + +If you are only interested in using RangeShifter, you can ignore this and head to the repo of one of the interfaces: + +- [WIP] RangeShifter GUI + +- [RangeShiftR](https://github.com/RangeShifter/RangeShiftR-pkg) + +- [RangeShifter-batch](https://github.com/RangeShifter/RangeShifter_batch) + +## Usage: git subtree + +In order to ensure that the same version of RangeShifter's core code is used by all three interfaces (RangeShiftR, RangeShifter-batch and the GUI), each interface repo keeps a copy of RScore as a git subtree. In this section, we describe how to use the git subtrees to update the subfolder and copy this repo anew. + +First, in a local clone of one of the interface repos, add a remote named `RScore` pointing to the RScore repo. This will be convenient as a shortcut for git subtree commands. + +```bash +git remote add RScore https://github.com/RangeShifter/RScore.git +``` + +### Pulling new changes + +To update the RScore subfolder with new changes made to the RScore repo, one can use the `git subtree pull` command: + +```bash +git subtree pull --prefix RScore +``` + +Note the path must match the location of the RScore subfolder, and the branch must match the one the subtree was originally added from (by default, this should be `main`). + +e.g. for RangeShifter-batch, use: + +```bash +git subtree pull --prefix src/RScore RScore main +``` + +while for RangeShiftR, use: + +```bash +git subtree pull --prefix RangeShiftR/src/RScore RScore main +``` + +### Pushing new changes to RScore + +We haven't yet found a way to push new changes made in a RScore subfolder back into the RScore repo. This is why we ask that contributions are made directly inside the RScore repo. + +If you know how to do this, please drop us a line! + +Alternatively, if you have already made changes on the subfolder, you could copy its contents into a new branch in RScore, then open a pull request. + +### Switching the subfolder to a new branch + +There is unfortunately to do so. To track a different branch of RScore, one must delete the RScore subfolder (via git) and import the subtree again: + +```bash +git rm src/RScore -r +git commit -m "switching subtree branch" +git subtree add --prefix src/RScore RScore +``` + +## Contributing to RangeShifter core code + +See [CONTRIBUTING](https://github.com/RangeShifter/RScore/blob/main/CONTRIBUTING.md). + +## Maintainers + +- [@JetteReeg](https://github.com/JetteReeg) +- [@TheoPannetier](https://github.com/TheoPannetier) diff --git a/src/RScore/RS_repos.png b/src/RScore/RS_repos.png new file mode 100644 index 0000000..b138a01 Binary files /dev/null and b/src/RScore/RS_repos.png differ diff --git a/src/RScore/RScore_logo.png b/src/RScore/RScore_logo.png new file mode 100644 index 0000000..dfce267 Binary files /dev/null and b/src/RScore/RScore_logo.png differ diff --git a/src/RScore/RSrandom.cpp b/src/RScore/RSrandom.cpp index d7153b0..ed50b31 100644 --- a/src/RScore/RSrandom.cpp +++ b/src/RScore/RSrandom.cpp @@ -26,7 +26,6 @@ #if !RS_RCPP - #if RSDEBUG #include "Parameters.h" extern paramSim* paramsSim; @@ -221,9 +220,11 @@ float RSrandom::FRandom(float min, float max) { return unif(*gen); } -int RSrandom::Bernoulli(double p) { - return Random() < p; -} + int RSrandom::Bernoulli(double p) { + if (p < 0) throw runtime_error("Bernoulli's p cannot be negative.\n"); + if (p > 1) throw runtime_error("Bernoulli's p cannot be above 1.\n"); + return Random() < p; + } double RSrandom::Normal(double mean, double sd) { return mean + sd * pNormal->operator()(*gen); @@ -303,4 +304,28 @@ double RSrandom::Cauchy(double loc, double scale) { #endif // RS_RCPP + +#if RSDEBUG + void testRSrandom() { + + { + // Bernoulli distribution + // Abuse cases + assert_error("Bernoulli's p cannot be negative.\n", []{ + RSrandom rsr; + rsr.Bernoulli(-0.3); + }); + assert_error("Bernoulli's p cannot be above 1.\n", [] { + RSrandom rsr; + rsr.Bernoulli(1.1); + }); + // Use cases + RSrandom rsr; + assert(rsr.Bernoulli(0) == 0); + assert(rsr.Bernoulli(1) == 1); + int bern_trial = rsr.Bernoulli(0.5); + assert(bern_trial == 0 || bern_trial == 1); + } + } +#endif // RSDEBUG //--------------------------------------------------------------------------- diff --git a/src/RScore/RSrandom.h b/src/RScore/RSrandom.h index 64e718e..6349520 100644 --- a/src/RScore/RSrandom.h +++ b/src/RScore/RSrandom.h @@ -38,6 +38,8 @@ #include #include +#include +#include "Utils.h" using namespace std; @@ -45,10 +47,7 @@ using namespace std; extern ofstream DEBUGLOG; #endif - - #if !RS_RCPP - //--------------- 2.) New version of RSrandom.cpp @@ -132,6 +131,9 @@ class RSrandom { #endif // !RS_RCPP +#if RSDEBUG + void testRSrandom(); +#endif // RSDEBUG //--------------------------------------------------------------------------- diff --git a/src/RScore/SNPTrait.cpp b/src/RScore/SNPTrait.cpp index ebda79a..102d9cf 100644 --- a/src/RScore/SNPTrait.cpp +++ b/src/RScore/SNPTrait.cpp @@ -21,14 +21,14 @@ SNPTrait::SNPTrait(SpeciesTrait* P) if (mutationDistribution != SSM && mutationDistribution != KAM) cout << endl << ("Error:: wrong mutation distribution for neutral markers, must be KAM or SSM \n"); - if (!mutationParameters.count(MAX)) + if (mutationParameters.count(MAX) != 1) cout << endl << ("Error:: KAM or SSM mutation distribution parameter must contain max value (e.g. max= ), max cannot exceed 256 \n"); if (wildType == -999) wildType = (int)mutationParameters.find(MAX)->second - 1; if (wildType > maxSNPAlleles - 1) - cout << endl << ("Error:: max number of alleles cannot exceed " << maxSNPAlleles << ".\n"); + cout << endl << "Error:: max number of alleles cannot exceed " << maxSNPAlleles << ".\n"; DistributionType initialDistribution = pProtoTrait->getInitialDistribution(); map initialParameters = pProtoTrait->getInitialParameters(); @@ -39,12 +39,12 @@ SNPTrait::SNPTrait(SpeciesTrait* P) switch (initialDistribution) { case UNIFORM: { - if (!initialParameters.count(MAX)) - cout << endl << ("Error:: initial SNP/Microsat distribution parameter must contain max value if set to UNIFORM (e.g. max= ), max cannot exceed " << maxSNPAlleles << "\n"); + if (initialParameters.count(MAX) != 1) + cout << endl << "Error:: initial SNP/Microsat distribution parameter must contain one max value if set to UNIFORM (e.g. max= ), max cannot exceed " << maxSNPAlleles << "\n"; float maxD = initialParameters.find(MAX)->second; if (maxD > maxSNPAlleles) { - cout << endl << ("Warning:: initial SNP/Microsat distribution parameter max cannot exceed " << maxSNPAlleles << ", resetting to " << maxSNPAlleles << "\n"); + cout << endl << "Warning:: initial SNP/Microsat distribution parameter max cannot exceed " << maxSNPAlleles << ", resetting to " << maxSNPAlleles << "\n"; maxD = maxSNPAlleles; //reserve 255 for wildtype } @@ -83,7 +83,6 @@ SNPTrait::SNPTrait(const SNPTrait& T) : // ---------------------------------------------------------------------------------------- void SNPTrait::mutate_KAM() { - const int positionsSize = pProtoTrait->getPositionsSize(); const auto& positions = pProtoTrait->getPositions(); const short ploidy = pProtoTrait->getPloidy(); @@ -199,7 +198,7 @@ void SNPTrait::inheritDiploid(sex_t whichChromosome, map nextBreakpoint) { std::advance(it, 1); nextBreakpoint = *it; - parentChromosome = !parentChromosome; //switch chromosome + parentChromosome = 1 - parentChromosome; //switch chromosome } if (locus <= nextBreakpoint) { diff --git a/src/RScore/SNPTrait.h b/src/RScore/SNPTrait.h index d4dd467..2133dbd 100644 --- a/src/RScore/SNPTrait.h +++ b/src/RScore/SNPTrait.h @@ -15,7 +15,7 @@ class SNPTrait : public TTrait { private: inline static int wildType = -999; - constexpr int maxSNPAlleles = UCHAR_MAX + 1; // allele is char, can take value 0-255 + const int maxSNPAlleles = UCHAR_MAX + 1; // allele is char, can take value 0-255 SpeciesTrait* pProtoTrait; diff --git a/src/RScore/Species.cpp b/src/RScore/Species.cpp index 81dd845..459c3d3 100644 --- a/src/RScore/Species.cpp +++ b/src/RScore/Species.cpp @@ -57,7 +57,6 @@ Species::Species(void) d0[i][j] = 0.0; alphaEmig[i][j] = 0.0; betaEmig[i][j] = 1.0; } } - // initialise transfer parameters moveModel = false; stgDepTrfr = false; sexDepTrfr = false; distMort = false; indVarTrfr = false; @@ -79,7 +78,6 @@ Species::Species(void) habDimTrfr = 0; straigtenPath = false; fullKernel = false; - // initialise settlement parameters stgDepSett = false; sexDepSett = false; indVarSett = false; minSteps = 0; maxSteps = 99999999; @@ -90,7 +88,6 @@ Species::Species(void) s0[i][j] = 1.0; alphaS[i][j] = 0.0; betaS[i][j] = 1.0; } } - // initialise attribute defaults spNum = 0; resetGeneticParameters(); @@ -114,7 +111,7 @@ short Species::getSpNum(void) { return spNum; } void Species::setDemogr(const demogrParams d) { if (d.repType >= 0 && d.repType <= 2) repType = d.repType; - if (d.repType >= 1) isDiploid = true; + if (d.repType >= 1) diploid = true; if (d.repSeasons >= 1) repSeasons = d.repSeasons; stageStruct = d.stageStruct; if (d.propMales > 0.0 && d.propMales < 1.0) propMales = d.propMales; @@ -360,8 +357,6 @@ void Species::deleteDDwtSurv(void) { delete[] ddwtSurv; ddwtSurv = 0; } } - -// Functions to handle min/max R or K (under environmental stochasticity) //void Species::setMinMax(float min,float max) { void Species::setMinMax(float min, float max) { if (min >= 0.0 && max > min) { @@ -375,8 +370,6 @@ float Species::getMinMax(short opt) { //--------------------------------------------------------------------------- -// Genetic functions - void Species::turnOffMutations(void) { mutationsOn = false; } @@ -395,12 +388,10 @@ void Species::resetGeneticParameters() { nIndsToSample = -9999; chromosomeEnds.clear(); samplePatchList.clear(); - stagesToSampleFrom.clear(); - } bool Species::isDiploid() const { - return isDiploid; + return diploid; } void Species::setNumberOfNeutralLoci(int nN) @@ -487,7 +478,7 @@ set Species::getTraitTypes() { } int Species::getNTraits() const { - return spTraitTable.size(); + return static_cast(spTraitTable.size()); } int Species::getNPositionsForTrait(const TraitType trait) const { @@ -506,19 +497,14 @@ set Species::getChromosomeEnds() const { return chromosomeEnds; } - //--------------------------------------------------------------------------- // Emigration functions - void Species::setEmig(const emigRules e) { -#if RSDEBUG - //DebugGUI("Species::setEmig(): e.indVar=" + Int2Str((int)e.indVar)); -#endif densDepEmig = e.densDep; stgDepEmig = e.stgDep; sexDepEmig = e.sexDep; indVarEmig = e.indVar; if (e.emigStage >= 0) emigStage = e.emigStage; - //setGenome(); +//setGenome(); } emigRules Species::getEmig(void) { @@ -555,8 +541,6 @@ float Species::getEmigD0(short stg, short sex) { } } - - //--------------------------------------------------------------------------- // Transfer functions @@ -570,7 +554,7 @@ void Species::setTrfr(const trfrRules t) { twinKern = t.twinKern; habMort = t.habMort; moveType = t.moveType; costMap = t.costMap; - //setGenome(); +//setGenome(); } trfrRules Species::getTrfr(void) { @@ -772,14 +756,10 @@ void Species::setSettTraits(const short stg, const short sex, const settleTraits } } -void Species::setGeneticParameters(const set& chromosomeEnds, const int genomeSize, const float recombinationRate, - const set& samplePatchList, const string nIndsToSample, const set& stagesToSampleFrom, string nSampleCellsFst) +void Species::setGeneticParameters(const std::set& chromosomeEnds, const int genomeSize, const float recombinationRate, + const std::set& samplePatchList, const string nIndsToSample, const std::set& stagesToSampleFrom, string nSampleCellsFst) { this->genomeSize = genomeSize; - for (auto position : chromosomeEnds) { - if (position > this->getGenomeSize() - 1) - cout << endl << "Genetics file: ERROR - chromosome ends " << position << " must not exceed genome size" << endl; - } this->chromosomeEnds = chromosomeEnds; this->recombinationRate = recombinationRate; this->samplePatchList = samplePatchList; @@ -802,9 +782,6 @@ settleTraits Species::getSettTraits(short stg, short sex) { return dd; } - - - //--------------------------------------------------------------------------- //--------------------------------------------------------------------------- //--------------------------------------------------------------------------- diff --git a/src/RScore/Species.h b/src/RScore/Species.h index 4666ae7..bedfcb7 100644 --- a/src/RScore/Species.h +++ b/src/RScore/Species.h @@ -31,9 +31,9 @@ The class holds all the demographic and dispersal parameters of the species. 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 @@ -49,6 +49,10 @@ #include "SpeciesTrait.h" #include "TTrait.h" #include +#include +#include + +class SpeciesTrait; // structures for demographic parameters @@ -280,7 +284,7 @@ class Species { short // option: 0 = return minimum, otherwise = return maximum ); - set& getSamplePatches() { + std::set& getSamplePatches() { return samplePatchList; }; @@ -288,7 +292,7 @@ class Species { return nIndsToSample; }; - set& getStagesToSample() { + std::set& getStagesToSample() { return stagesToSampleFrom; } @@ -419,19 +423,16 @@ class Species { //map>& getTraitTable(void); //return by reference so ensure variable recieving is const - set getTraitTypes(); + std::set getTraitTypes(); int getNTraits() const; int getNPositionsForTrait(const TraitType trait) const; int getGenomeSize() const; - void setGenomeSize(int); float getRecombinationRate() const; - void setRecombinationRate(float); - set getChromosomeEnds() const; - void setChromosomeEnds(const set& ends); - void setGeneticParameters(const set& chromosomeEnds, const int genomeSize, const float recombinationRate, - const set& samplePatchList, const string nIndsToSample, const set& stagesToSampleFrom, string nSampleCellsFst); - void setSamplePatchList(const set& samplePatchList); + std::set getChromosomeEnds() const; + void setGeneticParameters(const std::set& chromosomeEnds, const int genomeSize, const float recombinationRate, + const std::set& samplePatchList, const string nIndsToSample, const std::set& stagesToSampleFrom, string nSampleCellsFst); + void setSamplePatchList(const std::set& samplePatchList); private: @@ -484,17 +485,17 @@ class Species { // genome parameters /**The traits table.*/ - map> spTraitTable; - set chromosomeEnds; + std::map> spTraitTable; + std::set chromosomeEnds; int genomeSize; - bool isDiploid; + bool diploid; bool mutationsOn; int numberOfNeutralLoci; int numberOfAdaptiveTraits; float recombinationRate; - set samplePatchList; + std::set samplePatchList; string nSampleCellsFst; //for cell based landscape - set stagesToSampleFrom; + std::set stagesToSampleFrom; string nIndsToSample; //could be integer or 'all', all means in in selected patches not necessarily all in population // emigration parameters @@ -569,7 +570,6 @@ class Species { float betaS[NSTAGES][NSEXES]; // inflection point of the settlement reaction norm to density // other attributes - int spNum; }; diff --git a/src/RScore/SpeciesTrait.cpp b/src/RScore/SpeciesTrait.cpp index 7d5264a..f8c5d8a 100644 --- a/src/RScore/SpeciesTrait.cpp +++ b/src/RScore/SpeciesTrait.cpp @@ -1,6 +1,5 @@ #include "SpeciesTrait.h" -#include "Species.h" //could be handled in header file but here for now for flexibility SpeciesTrait::SpeciesTrait(vector parameters, Species* pSpecies) { @@ -20,7 +19,7 @@ SpeciesTrait::SpeciesTrait(vector parameters, Species* pSpecies) { if (traitType == SNP || traitType == ADAPTIVE) this->inherited = true; else - this->inherited = iequals(parameters[10], "true") ? true : false; + this->inherited = (parameters[10] == "true") ? true : false; if (this->isInherited()) { this->mutationDistribution = stringToDistributionType(parameters[11]); @@ -35,7 +34,7 @@ SpeciesTrait::SpeciesTrait(vector parameters, Species* pSpecies) { if (pSpecies->getNumberOfNeutralLoci() > 0) cout << endl << "Traits file: WARNING - can only have one set of neutral markers, overwriting previous" << endl; - else pSpecies->setNumberOfNeutralLoci(positions.size()); + else pSpecies->setNumberOfNeutralLoci(static_cast(positions.size())); } } @@ -53,6 +52,7 @@ TraitType SpeciesTrait::stringToTraitType(const std::string& str, sex_t sex) con else if (str == "kernel_probability") return KERNEL_PROBABILITY_M; else if (str == "crw_stepLength") return CRW_STEPLENGTH_M; else if (str == "crw_stepCorrelation") return CRW_STEPCORRELATION_M; + else throw logic_error(str + " is not a valid trait type."); } else { if (str == "emigration_d0") return E_D0_F; else if (str == "emigration_alpha") return E_ALPHA_F; @@ -71,6 +71,7 @@ TraitType SpeciesTrait::stringToTraitType(const std::string& str, sex_t sex) con else if (str == "sms_betaDB") return SMS_BETADB; else if (str == "neutral") return SNP; else if (str == "adaptive") return ADAPTIVE; + else throw logic_error(str + " is not a valid trait type."); } } @@ -79,6 +80,7 @@ ExpressionType SpeciesTrait::stringToExpressionType(const std::string& str) cons else if (str == "additive") return ADDITIVE; else if (str == "multiplicative") return MULTIPLICATIVE; else if (str == "#") return NEUTRAL; + else throw logic_error(str + " is not a valid gene expression type."); } DistributionType SpeciesTrait::stringToDistributionType(const std::string& str) const @@ -91,6 +93,7 @@ DistributionType SpeciesTrait::stringToDistributionType(const std::string& str) else if (str == "negExp") return NEGEXP; else if (str == "KAM") return KAM; else if (str == "SSM") return SSM; + else throw logic_error(str + " is not a valid distribution type."); } map SpeciesTrait::stringToParameterMap(string parameters) const { @@ -135,29 +138,38 @@ set SpeciesTrait::stringToLoci(string pos, string nLoci, Species* pSpecies) set positions; if (pos != "random") { - //stringstream ss(pos); - - //string value, valueWithin; - //while (std::getline(ss, value, ',')) { - // stringstream sss(value); - // vector minMax; - // while (std::getline(sss, valueWithin, '-')) { - // minMax.push_back(stoi(valueWithin)); - // } - // if (minMax[0] > pSpecies->getGenomeSize() || minMax[1] > pSpecies->getGenomeSize()) { - // cout << endl << "Traits file: ERROR - trait positions must not exceed genome size" << endl; - // } - // else { - // if (minMax.size() > 1) { - // for (int i = minMax[0]; i < minMax[1] + 1; ++i) { - // positions.insert(i); - // } - // } - // else - // positions.insert(minMax[0]); - // } - //} - positions = convertStringToSet(pos); + + // Parse comma-separated list from input string + stringstream ss(pos); + string value, valueWithin; + // Read comma-separated positions + while (std::getline(ss, value, ',')) { + stringstream sss(value); + vector positionRange; + // Read single positions and dash-separated ranges + while (std::getline(sss, valueWithin, '-')) { + positionRange.push_back(stoi(valueWithin)); + } + switch (positionRange.size()) + { + case 1: // single position + if (positionRange[0] > pSpecies->getGenomeSize()) + throw logic_error("Traits file: ERROR - trait positions must not exceed genome size"); + positions.insert(positionRange[0]); + break; + case 2: // dash-separated range + if (positionRange[0] > pSpecies->getGenomeSize() || positionRange[1] > pSpecies->getGenomeSize()) { + throw logic_error("Traits file: ERROR - trait positions must not exceed genome size"); + } + for (int i = positionRange[0]; i < positionRange[1] + 1; ++i) { + positions.insert(i); + } + break; + default: // zero or more than 2 values between commas: error + throw logic_error("Traits file: ERROR - incorrectly formatted position range."); + break; + } + } for (auto position : positions) { if (position > pSpecies->getGenomeSize()) diff --git a/src/RScore/SpeciesTrait.h b/src/RScore/SpeciesTrait.h index 3894207..9a634cb 100644 --- a/src/RScore/SpeciesTrait.h +++ b/src/RScore/SpeciesTrait.h @@ -2,6 +2,8 @@ #define SPECIESTRAITH #include "Parameters.h" +#include "Species.h" + #include #include #include @@ -34,7 +36,7 @@ class SpeciesTrait { float getMutationRate() const { return mutationRate; } short getPloidy() const { return ploidy; } set& getPositions() { return positions; } // returning by reference, make sure receiver is const - int getPositionsSize() const { return positions.size(); } + int getPositionsSize() const { return static_cast(positions.size()); } bool isInherited() const { return inherited; } DistributionType getMutationDistribution() const { return mutationDistribution; }; map getMutationParameters() const { return mutationParameters; }; diff --git a/src/RScore/SubCommunity.cpp b/src/RScore/SubCommunity.cpp index f0deddd..5164444 100644 --- a/src/RScore/SubCommunity.cpp +++ b/src/RScore/SubCommunity.cpp @@ -149,12 +149,6 @@ void SubCommunity::initialInd(Landscape* pLandscape, Species* pSpecies, else { if (iind.sex == 1) probmale = 1.0; else probmale = 0.0; } - //if (trfr.moveModel) { - // movt = true; - // if (trfr.moveType) { - // - // } - //} pInd = new Individual(pCell, pPatch, stg, age, repInt, probmale, trfr.moveModel, trfr.moveType); // add new individual to the population @@ -233,10 +227,6 @@ void SubCommunity::resetPopns(void) { void SubCommunity::resetPossSettlers(void) { if (subCommNum == 0) return; // not applicable in the matrix - //int npops = (int)popns.size(); - //for (int i = 0; i < npops; i++) { // all populations - // popns[i]->resetPossSettlers(); - //} pPatch->resetPossSettlers(); } @@ -294,7 +284,6 @@ void SubCommunity::reproduction(int resol, float epsGlobal, short rasterType, bo { if (subCommNum == 0) return; // no reproduction in the matrix float localK, envval; - //Species *pSpecies; Cell* pCell; envGradParams grad = paramsGrad->getGradient(); envStochParams env = paramsStoch->getStoch(); @@ -350,21 +339,6 @@ void SubCommunity::reproduction(int resol, float epsGlobal, short rasterType, bo } } */ - - /* - #if RSDEBUG - //if (npops > 0) { - // DEBUGLOG << "SubCommunity::reproduction(): this = " << this - // << " npops = " << npops << endl; - //} - #endif - #if RSDEBUG - //DEBUGLOG << "SubCommunity::reproduction(): patchNum = " << pPatch->getPatchNum() - // << " nCells " << pPatch->getNCells() - // << " localK = " << localK - // << endl; - #endif - */ } void SubCommunity::emigration(void) @@ -378,7 +352,6 @@ void SubCommunity::emigration(void) // NOTE that even if K is zero, it could have been >0 in previous time-step, and there // might be emigrants if there is non-juvenile emigration for (int i = 0; i < npops; i++) { // all populations - // localK = pPatch->getK(); popns[i]->emigration(localK); } } @@ -386,28 +359,13 @@ void SubCommunity::emigration(void) // Remove emigrants from their natal patch and add to patch 0 (matrix) void SubCommunity::initiateDispersal(SubCommunity* matrix) { if (subCommNum == 0) return; // no dispersal initiation in the matrix -#if RSDEBUG - //DEBUGLOG << "SubCommunity::initiateDispersal(): this=" << this - // << " subCommNum=" << subCommNum - // << endl; -#endif popStats pop; disperser disp; int npops = (int)popns.size(); -#if RSDEBUG - //DEBUGLOG << "SubCommunity::initiateDispersal(): patchNum = " << patchNum - // << " npops " << npops - // << endl; -#endif for (int i = 0; i < npops; i++) { // all populations pop = popns[i]->getStats(); for (int j = 0; j < pop.nInds; j++) { -#if RSDEBUG - //DEBUGLOG << "SubCommunity::initiateDispersal(): i = " << i - // << " j " << j - // << endl; -#endif disp = popns[i]->extractDisperser(j); if (disp.yes) { // disperser - has already been removed from natal population // add to matrix population @@ -472,10 +430,6 @@ void SubCommunity::addTransferDataForInd(trfrData* avgTrfrData) { // Add an individual into the local population of its species in the patch void SubCommunity::recruit(Individual* pInd, Species* pSpecies) { -#if RSDEBUG - //DEBUGLOG << "SubCommunity::recruit(): this=" << this - // << endl; -#endif int npops = (int)popns.size(); for (int i = 0; i < npops; i++) { // all populations if (pSpecies == popns[i]->getSpecies()) { @@ -485,33 +439,21 @@ void SubCommunity::recruit(Individual* pInd, Species* pSpecies) { } // Transfer through the matrix - run for the matrix sub-community only -#if RS_RCPP // included also SEASONAL +#if RS_RCPP int SubCommunity::transfer(Landscape* pLandscape, short landIx, short nextseason) #else int SubCommunity::transfer(Landscape* pLandscape, short landIx) -#endif // SEASONAL || RS_RCPP +#endif // RS_RCPP { -#if RSDEBUG - //DEBUGLOG << "SubCommunity::transfer(): this=" << this - // << endl; -#endif int ndispersers = 0; int npops = (int)popns.size(); for (int i = 0; i < npops; i++) { // all populations -#if RS_RCPP // included also SEASONAL +#if RS_RCPP ndispersers += popns[i]->transfer(pLandscape, landIx, nextseason); #else ndispersers += popns[i]->transfer(pLandscape, landIx); -#endif // SEASONAL || RS_RCPP -#if RSDEBUG - //DEBUGLOG << "SubCommunity::transfer(): i = " << i - // << " this = " << this - // << " subCommNum = " << subCommNum - // << " npops = " << npops - // << " popns[i] = " << popns[i] - // << " ndispersers = " << ndispersers - // << endl; -#endif +#endif // RS_RCPP + } return ndispersers; } @@ -524,12 +466,6 @@ int SubCommunity::transfer(Landscape* pLandscape, short landIx) void SubCommunity::completeDispersal(Landscape* pLandscape, bool connect) { -#if RSDEBUG - //DEBUGLOG << "SubCommunity::completeDispersal(): this=" << this - // << endl; -#endif -//popStats pop; -//int popsize,subcomm; int popsize; disperser settler; Species* pSpecies; @@ -543,64 +479,20 @@ void SubCommunity::completeDispersal(Landscape* pLandscape, bool connect) for (int i = 0; i < npops; i++) { // all populations pSpecies = popns[i]->getSpecies(); popsize = popns[i]->getNInds(); -#if RSDEBUG - //DEBUGLOG << "SubCommunity::completeDispersal(): i=" << i - // << " subCommNum=" << subCommNum - // << " pPatch=" << pPatch << " patchNum=" << pPatch->getPatchNum() - // << " npops=" << npops - // << " popns[i]=" << popns[i] - // << " popsize=" << popsize - // << endl; -#endif for (int j = 0; j < popsize; j++) { bool settled; settler = popns[i]->extractSettler(j); settled = settler.yes; if (settled) { - // settler - has already been removed from matrix population - // in popns[i]->extractSettler() -#if RSDEBUG -//locn loc = pNewCell->getLocn(); -//DEBUGLOG << "SubCommunity::completeDispersal(): j=" << j << " settled=" << settled; -//#if GROUPDISP -//if (emig.groupdisp) { -//DEBUGLOG << " groupID=" << groupsettler.pGroup->getGroupID() -// << " groupsize=" << groupsettler.groupsize -// << " status=" << groupsettler.pGroup->getStatus(); -//} -//else { -//DEBUGLOG << " settler.pInd=" << settler.pInd; -//} -//#else -//DEBUGLOG << " settler.pInd=" << settler.pInd; -//#endif // GROUPDISP -//DEBUGLOG << " loc.x=" << loc.x << " loc.y=" << loc.y << endl; -#endif // RSDEBUG + // settler - has already been removed from matrix population // find new patch pNewPatch = (Patch*)settler.pCell->getPatch(); // find population within the patch (if there is one) - // subcomm = ppatch->getSubComm(); - // if (subcomm == 0) { // no sub-community has yet been set up - // // CANNOT SET UP A NEW ONE FROM WITHIN AN EXISTING ONE!!!!! - // } pPop = (Population*)pNewPatch->getPopn((intptr)pSpecies); -#if RSDEBUG - //DEBUGLOG << "SubCommunity::completeDispersal(): j = " << j - // << " pCell = " << pCell << " ppatch = " << ppatch << " pPop = " << pPop - // << endl; -#endif if (pPop == 0) { // settler is the first in a previously uninhabited patch // create a new population in the corresponding sub-community pSubComm = (SubCommunity*)pNewPatch->getSubComm(); -#if RSDEBUG - //DEBUGLOG << "SubCommunity::completeDispersal(): j = " << j - // << " pSubComm = " << pSubComm << endl; -#endif pPop = pSubComm->newPopn(pLandscape, pSpecies, pNewPatch, 0); -#if RSDEBUG - //DEBUGLOG << "SubCommunity::completeDispersal(): j=" << j - // << " pPop=" << pPop << endl; -#endif } pPop->recruit(settler.pInd); if (connect) { // increment connectivity totals @@ -619,16 +511,6 @@ void SubCommunity::completeDispersal(Landscape* pLandscape, bool connect) } // remove pointers in the matrix popn to settlers popns[i]->clean(); -#if RSDEBUG - //pop = popns[i]->getStats(); - popsize = popns[i]->getNInds(); - //DEBUGLOG << "SubCommunity::completeDispersal(): i=" << i - // << " popns[i]=" << popns[i] - //// << " pop.pPatch = " << pop.pPatch - //// << " pop.nInds = " << pop.nInds - // << " popsize=" << popsize - // << endl; -#endif } } @@ -692,10 +574,6 @@ void SubCommunity::createOccupancy(int nrows) { } void SubCommunity::updateOccupancy(int row) { -#if RSDEBUG - //DEBUGLOG << "SubCommunity::updateOccupancy(): this=" << this - // << endl; -#endif popStats pop; int npops = (int)popns.size(); for (int i = 0; i < npops; i++) { @@ -762,7 +640,6 @@ void SubCommunity::outPop(Landscape* pLandscape, int rep, int yr, int gen) // or it is the matrix patch in a patch-based model int npops = (int)popns.size(); int patchnum; - //Species* pSpecies; Cell* pCell; float localK; float eps = 0.0; @@ -778,7 +655,6 @@ void SubCommunity::outPop(Landscape* pLandscape, int rep, int yr, int gen) patchnum = pPatch->getPatchNum(); for (int i = 0; i < npops; i++) { // all populations - // pSpecies = popns[i]->getSpecies(); localK = pPatch->getK(); if (localK > 0.0 || (land.patchModel && patchnum == 0)) { popns[i]->outPopulation(rep, yr, gen, eps, land.patchModel, writeEnv, gradK); @@ -829,7 +705,6 @@ int SubCommunity::stagePop(int stage) { // Open traits file and write header record bool SubCommunity::outTraitsHeaders(Landscape* pLandscape, Species* pSpecies, int landNr) { - //Population *pPop; landParams land = pLandscape->getLandParams(); if (landNr == -999) { // close file if (outtraits.is_open()) outtraits.close(); @@ -1236,9 +1111,6 @@ traitsums SubCommunity::outTraits(traitCanvas tcanv, } // CURRENTLY INDIVIDUAL VARIATION CANNOT BE SEX-DEPENDENT -// ngenes = 1; - double mnS0[2], mnAlpha[2], mnBeta[2], sdS0[2], sdAlpha[2], sdBeta[2]; - if (writefile) outtraits << endl; for (int s = 0; s < NSEXES; s++) { diff --git a/src/RScore/SubCommunity.h b/src/RScore/SubCommunity.h index 2d211af..19017b0 100644 --- a/src/RScore/SubCommunity.h +++ b/src/RScore/SubCommunity.h @@ -32,9 +32,9 @@ CURRENTLY the number of Populations withn a SubCommunity is LIMITED TO ONE. 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 @@ -105,7 +105,7 @@ class SubCommunity { Individual*, // pointer to Individual Species* // pointer to Species ); -#if RS_RCPP // included also SEASONAL +#if RS_RCPP int transfer( // Transfer through matrix - run for matrix SubCommunity only Landscape*, // pointer to Landscape short, // landscape change index @@ -116,7 +116,7 @@ class SubCommunity { Landscape*, // pointer to Landscape short // landscape change index ); -#endif // SEASONAL || RS_RCPP +#endif // RS_RCPP // Remove emigrants from patch 0 (matrix) and transfer to SubCommunity in which // their destination co-ordinates fall (executed for the matrix patch only) void completeDispersal( @@ -194,10 +194,9 @@ class SubCommunity { private: intptr subCommNum; // SubCommunity number - // 0 is reserved for the SubCommunity in the inter-patch matrix -// intptr *occupancy; // pointer to occupancy array - Patch* pPatch; - int* occupancy; // pointer to occupancy array + // 0 is reserved for the SubCommunity in the inter-patch matrix + Patch *pPatch; + int *occupancy; // pointer to occupancy array std::vector popns; bool initial; // WILL NEED TO BE CHANGED FOR MULTIPLE SPECIES ... diff --git a/src/RScore/TraitFactory.h b/src/RScore/TraitFactory.h index 9cf1dda..e7420ab 100644 --- a/src/RScore/TraitFactory.h +++ b/src/RScore/TraitFactory.h @@ -4,7 +4,7 @@ #include "SpeciesTrait.h" #include "SNPTrait.h" #include "QTLTrait.h" -#include "AdaptiveTrait.h" +#include "GeneticLoad.h" class TraitFactory { diff --git a/src/RScore/Utils.cpp b/src/RScore/Utils.cpp new file mode 100644 index 0000000..9d36b92 --- /dev/null +++ b/src/RScore/Utils.cpp @@ -0,0 +1,26 @@ +#include "Utils.h" + +const string Int2Str(const int x) +{ + ostringstream o; + if (!(o << x)) return "ERROR"; + return o.str(); +} +const string Float2Str(const float x) { + ostringstream o; + if (!(o << x)) return "ERROR"; + return o.str(); +} +const string Double2Str(const double x) { + ostringstream o; + if (!(o << x)) return "ERROR"; + return o.str(); +} + +// Evaluate a lambda and assert we get the correct error +void assert_error(const string& exptd_err_msg, void (*lambda)(void)) { + string err_msg{ "No error.\n" }; + try { lambda(); } // evaluate + catch (exception& e) { err_msg = e.what(); } + assert(err_msg == exptd_err_msg); +} \ No newline at end of file diff --git a/src/RScore/Utils.h b/src/RScore/Utils.h new file mode 100644 index 0000000..34854a7 --- /dev/null +++ b/src/RScore/Utils.h @@ -0,0 +1,18 @@ +#ifndef UtilsH +#define UtilsH + +#include +#include +#include + +#include +using namespace std; + +const string Int2Str(const int x); +const string Float2Str(const float x); +const string Double2Str(const double x); + +// Evaluate a lambda and assert we get the correct error +void assert_error(const string& exptd_err_msg, void (*x)(void)); + +#endif // UtilsH diff --git a/src/RScore/branches.png b/src/RScore/branches.png new file mode 100644 index 0000000..12e3a96 Binary files /dev/null and b/src/RScore/branches.png differ diff --git a/src/RScore/git_cheatsheet.md b/src/RScore/git_cheatsheet.md new file mode 100644 index 0000000..53ca175 --- /dev/null +++ b/src/RScore/git_cheatsheet.md @@ -0,0 +1,107 @@ +# Git Cheatsheet + +Quick reference on Git usage for RangeShifter contributors + +#### Creating a local copy of this repo + +```bash +git clone https://github.com/RangeShifter/RScore.git +``` + +#### Enquire about the current state of changes + +``` +git status +``` +This will display whether the local branch is up-to-date with its GitHub counterpart, what files have been changed and if they are staged for commit or not. + +#### Updating the active branch with latest changes (pulling) + +```bash +git pull +``` + +#### Updating the active branch with another branch (merging) + +```bash +git merge +``` + +Merging may trigger a merge conflict is the same lines have been changed on both branches. See Solving Merge Conflict below. + +#### Staging changes to be committed +Changes to local files must first be added to the commit queue ("staged") before being committed. + +```bash +git add # single file +git add . # entire active folder +``` + +#### Creating a commit + +```bash +git commit -m "commit message" +``` + +A good commit message is concise, but descriptive. + +#### Uploading local commits to GitHub (pushing) + +```bash +git push +``` + +#### Switching to an existing branch + +```bash +git checkout +``` +Switching does not update the new active branch with GitHub automatically, so make sure to pull after switching! + +#### Creating a new branch from active branch + +```bash +git branch +``` +You will also need to set the corresponding branch on `origin` (GitHub) before you can push: + +```bash +git push --set-upstream origin +``` + +New branches can also be created on GitHub (drop-down button in the top-left corner of the main page). +New branches on GitHub are brought to the local copy with the next pull. + +#### Deleting a branch locally + +```bash +git branch -d +``` + +#### Instruct Git to not track some files + +Open `.gitignore` and add the path to the files or folders to exclude from the git history. +Check with `git status` that the files are indeed not tracked by git anymore. + +#### Solving a merge conflict +Merge conflicts can arise when multiple contributors simulatneously change the same line of code. +In such cases, git is incapable of deciding which version should be kept and asks for human input. +Git tells you which files are affected by the conflict. +Open each file and resolve **each** section that looks like this: + +``` +After opening the box, we can affirm that +<<<<<<<<<<< HEAD # delete this line +Shroedinger's cat is alive. # delete either this... +================ # delete this line +Shroedinger's cat is dead. # ... or this (or keep both, or none, or a different solution) +>>>>>>>>>>> SHA # delete this line +What an insightful result! +``` + +Ctrl+F "HEAD" is really helpful for finding the conflicts in large files. +When you are done, create a commit stating e.g. "solved merge conflict" and push. + +## Git subtrees + +See [Git subtrees](https://github.com/RangeShifter/RScore/tree/development-guidelines#usage-git-subtree) section in README.