shift from taus2 and MT19937-64 to PCG RNGs

Author: bhaller

CMakeLists.txt

@@ -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()

SLiM.xcodeproj/project.pbxproj

@@ -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>"; };
@@ -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>";

SLiM.xcodeproj/xcshareddata/xcschemes/SLiM.xcscheme

@@ -40,7 +40,7 @@
       </Testables>
    </TestAction>
    <LaunchAction
-      buildConfiguration = "Debug"
+      buildConfiguration = "Release"
       selectedDebuggerIdentifier = "Xcode.DebuggerFoundation.Debugger.LLDB"
       selectedLauncherIdentifier = "Xcode.DebuggerFoundation.Launcher.LLDB"
       launchStyle = "0"
@@ -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;"
@@ -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;"
@@ -98,7 +98,7 @@
          </CommandLineArgument>
          <CommandLineArgument
             argument = "-testSLiM"
-            isEnabled = "NO">
+            isEnabled = "YES">
          </CommandLineArgument>
          <CommandLineArgument
             argument = "-x"

VERSIONS

@@ -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):

core/chromosome.cpp

@@ -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);
 	}
@@ -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;
@@ -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;
@@ -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
 
@@ -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" –
@@ -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
@@ -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);
 	}
@@ -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
@@ -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);
 		}
@@ -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);
 		}
@@ -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);