Skip to content

shift to PCG RNGs #571

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
bhaller merged 1 commit into from
Nov 1, 2025
Merged

shift to PCG RNGs #571

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -90,7 +90,7 @@ if(CPPCHECK)
find_program(CPPCHECK_EXECUTABLE cppcheck)
if (CPPCHECK_EXECUTABLE)
message(STATUS "CPPCHECK is ${CPPCHECK}; building with cppcheck (for development)")
set(CPPCHECK_COMMAND "${CPPCHECK_EXECUTABLE}" "--enable=all;--suppress=missingIncludeSystem;--suppress=syntaxError;--suppress=unmatchedSuppression;--inline-suppr;--std=c++11;--quiet;--suppress=*:robin_hood.h;--suppress=*:lodepng.cpp;--suppress=checkersReport;--suppress=*:eidos_openmp.h;--suppress=unusedFunction")
set(CPPCHECK_COMMAND "${CPPCHECK_EXECUTABLE}" "--enable=all;--suppress=missingIncludeSystem;--suppress=syntaxError;--suppress=unmatchedSuppression;--inline-suppr;--std=c++11;--quiet;--suppress=*:robin_hood.h;--suppress=*:lodepng.cpp;--suppress=*:pcg_extras.hpp;--suppress=*:pcg_random.hpp;--suppress=checkersReport;--suppress=*:eidos_openmp.h;--suppress=unusedFunction")
message(STATUS "+++ cppcheck is at ${CPPCHECK_EXECUTABLE}")
message(STATUS "+++ CPPCHECK_COMMAND is ${CPPCHECK_COMMAND}")
else()
Expand Down
4 changes: 4 additions & 0 deletions SLiM.xcodeproj/project.pbxproj
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2189,6 +2189,8 @@
98D7D6642AB24CBC002AFE34 /* chisq.c */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.c; path = chisq.c; sourceTree = "<group>"; };
98D7EBEE28CE557C00DEAAC4 /* eidos_multi */ = {isa = PBXFileReference; explicitFileType = "compiled.mach-o.executable"; includeInIndex = 0; path = eidos_multi; sourceTree = BUILT_PRODUCTS_DIR; };
98D7ED2D28CE58FC00DEAAC4 /* slim_multi */ = {isa = PBXFileReference; explicitFileType = "compiled.mach-o.executable"; includeInIndex = 0; path = slim_multi; sourceTree = BUILT_PRODUCTS_DIR; };
98D957882EB53494008314C1 /* pcg_extras.hpp */ = {isa = PBXFileReference; lastKnownFileType = sourcecode.cpp.h; path = pcg_extras.hpp; sourceTree = "<group>"; };
98D957892EB53494008314C1 /* pcg_random.hpp */ = {isa = PBXFileReference; lastKnownFileType = sourcecode.cpp.h; path = pcg_random.hpp; sourceTree = "<group>"; };
98DB3D6D1E6122AE00E2C200 /* interaction_type.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = interaction_type.cpp; sourceTree = "<group>"; };
98DB3D6E1E6122AE00E2C200 /* interaction_type.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = interaction_type.h; sourceTree = "<group>"; };
98DC9838289986B300160DD8 /* GitSHA1.cpp.in */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = text; path = GitSHA1.cpp.in; sourceTree = "<group>"; };
Expand Down Expand Up @@ -2389,6 +2391,8 @@
98235687252FDD120096A745 /* lodepng.h */,
98235681252FDCF50096A745 /* lodepng.cpp */,
98186DB8254A8B1600F9118C /* robin_hood.h */,
98D957882EB53494008314C1 /* pcg_extras.hpp */,
98D957892EB53494008314C1 /* pcg_random.hpp */,
);
name = dependencies;
sourceTree = "<group>";
Expand Down
8 changes: 4 additions & 4 deletions SLiM.xcodeproj/xcshareddata/xcschemes/SLiM.xcscheme
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,7 @@
</Testables>
</TestAction>
<LaunchAction
buildConfiguration = "Debug"
buildConfiguration = "Release"
selectedDebuggerIdentifier = "Xcode.DebuggerFoundation.Debugger.LLDB"
selectedLauncherIdentifier = "Xcode.DebuggerFoundation.Launcher.LLDB"
launchStyle = "0"
Expand Down Expand Up @@ -70,7 +70,7 @@
</CommandLineArgument>
<CommandLineArgument
argument = "-time"
isEnabled = "YES">
isEnabled = "NO">
</CommandLineArgument>
<CommandLineArgument
argument = "-d MUTDEFS=\&apos;initializeMutationType(\\&quot;m1\\&quot;, 0.5, &quot;f&quot;, 0.0);\&apos; &quot;/Users/bhaller/Desktop/cmdline test.slim&quot;"
Expand All @@ -86,7 +86,7 @@
</CommandLineArgument>
<CommandLineArgument
argument = "-mem"
isEnabled = "YES">
isEnabled = "NO">
</CommandLineArgument>
<CommandLineArgument
argument = "-define &quot;N=1000&quot; -d &quot;THETA=5&quot; -d &quot;mu=THETA/(4*N)&quot;"
Expand All @@ -98,7 +98,7 @@
</CommandLineArgument>
<CommandLineArgument
argument = "-testSLiM"
isEnabled = "NO">
isEnabled = "YES">
</CommandLineArgument>
<CommandLineArgument
argument = "-x"
Expand Down
1 change: 1 addition & 0 deletions VERSIONS
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@ development head (in the master branch):
extend text() to support drawing text and an angle, with new [float angle = 0.0] parameter
add mtext() call to Plot, for drawing text in the margins outside the plot area
add rowSums() and colSums() functions to Eidos, for use with matrices as a faster alternative to apply()
add the PCG random number generator, switch to pcg32_fast and pcg64_fast, remove all use of the old taus2 and MT19937-64 generators; note this completely breaks backward reproducibility


version 5.1 (Eidos version 4.1):
Expand Down
55 changes: 27 additions & 28 deletions core/chromosome.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -990,19 +990,17 @@ int Chromosome::DrawSortedUniquedMutationPositions(int p_count, IndividualSex p_
}

// draw all the positions, and keep track of the genomic element type for each
gsl_rng *rng = EIDOS_GSL_RNG(omp_get_thread_num());
Eidos_MT_State *mt = EIDOS_MT_RNG(omp_get_thread_num());
gsl_rng *rng_gsl = EIDOS_GSL_RNG(omp_get_thread_num());
EidosRNG_64_bit &rng_64 = EIDOS_64BIT_RNG(omp_get_thread_num());

for (int i = 0; i < p_count; ++i)
{
int mut_subrange_index = static_cast<int>(gsl_ran_discrete(rng, lookup));
int mut_subrange_index = static_cast<int>(gsl_ran_discrete(rng_gsl, lookup));
const GESubrange &subrange = (*subranges)[mut_subrange_index];
GenomicElement *source_element = subrange.genomic_element_ptr_;

// Draw the position along the chromosome for the mutation, within the genomic element
slim_position_t position = subrange.start_position_ + static_cast<slim_position_t>(Eidos_rng_uniform_int_MT64(mt, subrange.end_position_ - subrange.start_position_ + 1));
// old 32-bit position not MT64 code:
//slim_position_t position = subrange.start_position_ + static_cast<slim_position_t>(Eidos_rng_uniform_int(rng, (uint32_t)(subrange.end_position_ - subrange.start_position_ + 1)));
slim_position_t position = subrange.start_position_ + static_cast<slim_position_t>(Eidos_rng_interval_uint64(rng_64, subrange.end_position_ - subrange.start_position_ + 1));

p_positions.emplace_back(position, source_element);
}
Expand Down Expand Up @@ -1310,8 +1308,8 @@ MutationIndex Chromosome::DrawNewMutationExtended(std::pair<slim_position_t, Gen

// OK, now we know the background nucleotide; determine the mutation rates to derived nucleotides
double *nuc_thresholds = genomic_element_type.mm_thresholds + (size_t)original_nucleotide * 4;
gsl_rng *rng = EIDOS_GSL_RNG(omp_get_thread_num());
double draw = Eidos_rng_uniform(rng);
EidosRNG_64_bit &rng_64 = EIDOS_64BIT_RNG(omp_get_thread_num());
double draw = Eidos_rng_uniform_doubleCO(rng_64);

if (draw < nuc_thresholds[0]) nucleotide = 0;
else if (draw < nuc_thresholds[1]) nucleotide = 1;
Expand Down Expand Up @@ -1384,8 +1382,8 @@ MutationIndex Chromosome::DrawNewMutationExtended(std::pair<slim_position_t, Gen
// OK, now we know the background nucleotide; determine the mutation rates to derived nucleotides
int trinuc = ((int)background_nuc1) * 16 + ((int)original_nucleotide) * 4 + (int)background_nuc3;
double *nuc_thresholds = genomic_element_type.mm_thresholds + (size_t)trinuc * 4;
gsl_rng *rng = EIDOS_GSL_RNG(omp_get_thread_num());
double draw = Eidos_rng_uniform(rng);
EidosRNG_64_bit &rng_64 = EIDOS_64BIT_RNG(omp_get_thread_num());
double draw = Eidos_rng_uniform_doubleCO(rng_64);

if (draw < nuc_thresholds[0]) nucleotide = 0;
else if (draw < nuc_thresholds[1]) nucleotide = 1;
Expand Down Expand Up @@ -1486,13 +1484,13 @@ void Chromosome::_DrawCrossoverBreakpoints(IndividualSex p_parent_sex, const int
}

// draw recombination breakpoints
gsl_rng *rng = EIDOS_GSL_RNG(omp_get_thread_num());
Eidos_MT_State *mt = EIDOS_MT_RNG(omp_get_thread_num());
gsl_rng *rng_gsl = EIDOS_GSL_RNG(omp_get_thread_num());
EidosRNG_64_bit &rng_64 = EIDOS_64BIT_RNG(omp_get_thread_num());

for (int i = 0; i < p_num_breakpoints; i++)
{
slim_position_t breakpoint = 0;
int recombination_interval = static_cast<int>(gsl_ran_discrete(rng, lookup));
int recombination_interval = static_cast<int>(gsl_ran_discrete(rng_gsl, lookup));

// choose a breakpoint anywhere in the chosen recombination interval with equal probability

Expand All @@ -1508,7 +1506,7 @@ void Chromosome::_DrawCrossoverBreakpoints(IndividualSex p_parent_sex, const int
// positions to the left of its enclosed bases, up to and including the position to the left of the final base given as the
// end position of the interval. The next interval's first owned recombination position is therefore to the left of the
// base that is one position to the right of the end of the preceding interval. So we have to add one to the position
// given by recombination_end_positions_[recombination_interval - 1], at minimum. Since Eidos_rng_uniform_int() returns
// given by recombination_end_positions_[recombination_interval - 1], at minimum. Since Eidos_rng_interval_uint64() returns
// a zero-based random number, that means we need a +1 here as well.
//
// The key fact here is that a recombination breakpoint position of 1 means "break to the left of the base at position 1" –
Expand All @@ -1517,7 +1515,7 @@ void Chromosome::_DrawCrossoverBreakpoints(IndividualSex p_parent_sex, const int
// breakpoint. When their position is *equal*, the breakpoint gets serviced by switching strands. That logic causes the
// breakpoints to fall to the left of their designated base.
//
// Note that Eidos_rng_uniform_int() crashes (well, aborts fatally) if passed 0 for n. We need to guarantee that that doesn't
// Note that Eidos_rng_interval_uint64() crashes (well, aborts fatally) if passed 0 for n. We need to guarantee that that doesn't
// happen, and we don't want to waste time checking for that condition here. For a 1-base model, we are guaranteed that
// the overall recombination rate will be zero, by the logic in InitializeDraws(), and so we should not be called in the
// first place. For longer chromosomes that start with a 1-base recombination interval, the rate calculated by
Expand All @@ -1526,9 +1524,9 @@ void Chromosome::_DrawCrossoverBreakpoints(IndividualSex p_parent_sex, const int
// since we guarantee that recombination end positions are in strictly ascending order. So we should never crash. :->

if (recombination_interval == 0)
breakpoint = static_cast<slim_position_t>(Eidos_rng_uniform_int_MT64(mt, (*end_positions)[recombination_interval]) + 1);
breakpoint = static_cast<slim_position_t>(Eidos_rng_interval_uint64(rng_64, (*end_positions)[recombination_interval]) + 1);
else
breakpoint = (*end_positions)[recombination_interval - 1] + 1 + static_cast<slim_position_t>(Eidos_rng_uniform_int_MT64(mt, (*end_positions)[recombination_interval] - (*end_positions)[recombination_interval - 1]));
breakpoint = (*end_positions)[recombination_interval - 1] + 1 + static_cast<slim_position_t>(Eidos_rng_interval_uint64(rng_64, (*end_positions)[recombination_interval] - (*end_positions)[recombination_interval - 1]));

p_crossovers.emplace_back(breakpoint);
}
Expand Down Expand Up @@ -1604,7 +1602,9 @@ void Chromosome::_DrawDSBBreakpoints(IndividualSex p_parent_sex, const int p_num
// to a collision, because such redrawing would be liable to produce bias towards shorter extents. (Redrawing the crossover/
// noncrossover and simple/complex decisions would probably be harmless, but it is simpler to just make all decisions up front.)
int try_count = 0;
gsl_rng *rng = EIDOS_GSL_RNG(omp_get_thread_num());
gsl_rng *rng_gsl = EIDOS_GSL_RNG(omp_get_thread_num());
EidosRNG_64_bit &rng_64 = EIDOS_64BIT_RNG(omp_get_thread_num());

static std::vector<std::tuple<slim_position_t, slim_position_t, bool, bool>> dsb_infos; // using a static prevents reallocation

// If the redrawLengthsOnFailure parameter to initializeGeneConversion() is T, we jump back here on layout failure
Expand All @@ -1618,8 +1618,8 @@ void Chromosome::_DrawDSBBreakpoints(IndividualSex p_parent_sex, const int p_num
for (int i = 0; i < p_num_breakpoints; i++)
{
// If the gene conversion tract mean length is < 2.0, gsl_ran_geometric() will blow up, and we should treat the tract length as zero
bool noncrossover = (Eidos_rng_uniform(rng) <= non_crossover_fraction_); // tuple position 2
bool simple = (Eidos_rng_uniform(rng) <= simple_conversion_fraction_); // tuple position 3
bool noncrossover = (Eidos_rng_uniform_doubleCO(rng_64) <= non_crossover_fraction_); // tuple position 2
bool simple = (Eidos_rng_uniform_doubleCO(rng_64) <= simple_conversion_fraction_); // tuple position 3

dsb_infos.emplace_back(0, 0, noncrossover, simple);
}
Expand All @@ -1628,10 +1628,10 @@ void Chromosome::_DrawDSBBreakpoints(IndividualSex p_parent_sex, const int p_num
{
for (int i = 0; i < p_num_breakpoints; i++)
{
slim_position_t extent1 = gsl_ran_geometric(rng, gene_conversion_inv_half_length_); // tuple position 0
slim_position_t extent2 = gsl_ran_geometric(rng, gene_conversion_inv_half_length_); // tuple position 1
bool noncrossover = (Eidos_rng_uniform(rng) <= non_crossover_fraction_); // tuple position 2
bool simple = (Eidos_rng_uniform(rng) <= simple_conversion_fraction_); // tuple position 3
slim_position_t extent1 = gsl_ran_geometric(rng_gsl, gene_conversion_inv_half_length_); // tuple position 0
slim_position_t extent2 = gsl_ran_geometric(rng_gsl, gene_conversion_inv_half_length_); // tuple position 1
bool noncrossover = (Eidos_rng_uniform_doubleCO(rng_64) <= non_crossover_fraction_); // tuple position 2
bool simple = (Eidos_rng_uniform_doubleCO(rng_64) <= simple_conversion_fraction_); // tuple position 3

dsb_infos.emplace_back(extent1, extent2, noncrossover, simple);
}
Expand All @@ -1650,19 +1650,18 @@ void Chromosome::_DrawDSBBreakpoints(IndividualSex p_parent_sex, const int p_num
}

// First draw DSB points; dsb_points contains positions and a flag for whether the breakpoint is at a rate=0.5 position
Eidos_MT_State *mt = EIDOS_MT_RNG(omp_get_thread_num());
static std::vector<std::pair<slim_position_t, bool>> dsb_points; // using a static prevents reallocation
dsb_points.resize(0);

for (int i = 0; i < p_num_breakpoints; i++)
{
slim_position_t breakpoint = 0;
int recombination_interval = static_cast<int>(gsl_ran_discrete(rng, lookup));
int recombination_interval = static_cast<int>(gsl_ran_discrete(rng_gsl, lookup));

if (recombination_interval == 0)
breakpoint = static_cast<slim_position_t>(Eidos_rng_uniform_int_MT64(mt, (*end_positions)[recombination_interval]) + 1);
breakpoint = static_cast<slim_position_t>(Eidos_rng_interval_uint64(rng_64, (*end_positions)[recombination_interval]) + 1);
else
breakpoint = (*end_positions)[recombination_interval - 1] + 1 + static_cast<slim_position_t>(Eidos_rng_uniform_int_MT64(mt, (*end_positions)[recombination_interval] - (*end_positions)[recombination_interval - 1]));
breakpoint = (*end_positions)[recombination_interval - 1] + 1 + static_cast<slim_position_t>(Eidos_rng_interval_uint64(rng_64, (*end_positions)[recombination_interval] - (*end_positions)[recombination_interval - 1]));

if ((*rates)[recombination_interval] == 0.5)
dsb_points.emplace_back(breakpoint, true);
Expand Down
Loading
Loading