add a bunch of multitrait debug test code

Author: bhaller

core/individual.cpp

@@ -417,6 +417,13 @@ class Individual : public EidosDictionaryUnretained
 	template <const bool f_additiveTrait, const bool f_callbacks, const bool f_singlecallback>
 	void _IncorporateEffects_Diploid(Species *species, Haplosome *haplosome1, Haplosome *haplosome2, slim_trait_index_t trait_index, std::vector<SLiMEidosBlock*> &p_mutationEffect_callbacks);
 
+	// Debugging checkback for phenotype calculation; this is very slow, and does not use the non-neutral cache
+	slim_effect_t _CheckPhenotypeForTrait(slim_trait_index_t trait_index);
+	
+	void _Check_IncorporateEffects_Haploid(Species *species, Haplosome *haplosome, Trait *trait, std::vector<SLiMEidosBlock*> &p_mutationEffect_callbacks);
+	void _Check_IncorporateEffects_Hemizygous(Species *species, Haplosome *haplosome, Trait *trait, std::vector<SLiMEidosBlock*> &p_mutationEffect_callbacks);
+	void _Check_IncorporateEffects_Diploid(Species *species, Haplosome *haplosome1, Haplosome *haplosome2, Trait *trait, std::vector<SLiMEidosBlock*> &p_mutationEffect_callbacks);
+	
 	// for Subpopulation::ExecuteMethod_takeMigrants()
 	friend Subpopulation;
 };

core/individual.h

@@ -653,6 +653,30 @@ std::string StringForSLiMCycleStage(SLiMCycleStage p_stage)
 }
 
 // stream output for enumerations
+std::string StringForTraitType(TraitType p_trait_type)
+{
+	switch (p_trait_type)
+	{
+		case TraitType::kAdditive:								return gStr_additive;
+		case TraitType::kMultiplicative:						return gStr_multiplicative;
+	}
+	EIDOS_TERMINATION << "ERROR (StringForTraitType): (internal error) unexpected p_trait_type value." << EidosTerminate();
+}
+
+TraitType TraitTypeForString(std::string type)
+{
+	if (type == gStr_additive)	return TraitType::kAdditive;
+	else if (type == gStr_multiplicative)	return TraitType::kMultiplicative;
+	else
+		EIDOS_TERMINATION << "ERROR (TraitTypeForString): unrecognized triat type '" << type << "'." << EidosTerminate();
+}
+
+std::ostream& operator<<(std::ostream& p_out, TraitType p_trait_type)
+{
+	p_out << StringForTraitType(p_trait_type);
+	return p_out;
+}
+
 std::string StringForChromosomeType(ChromosomeType p_chromosome_type)
 {
 	switch (p_chromosome_type)
@@ -1629,6 +1653,8 @@ const std::string &gStr_willAutolog = EidosRegisteredString("willAutolog", gID_w
 const std::string &gStr_context = EidosRegisteredString("context", gID_context);
 
 // mostly other fixed strings
+const std::string gStr_additive = "additive"; // these trait type strings are not registered, no need
+const std::string gStr_multiplicative = "multiplicative";
 const std::string gStr_A = "A";	// these nucleotide strings are not registered, no need
 const std::string gStr_C = "C";
 const std::string gStr_G = "G";

core/slim_globals.cpp

@@ -492,6 +492,7 @@ void AccumulateMemoryUsageIntoTotal_Community(SLiMMemoryUsage_Community &p_usage
 #define DEBUG_BLOCK_REG_DEREG		0		// turn on to get logging about script block registration/deregistration
 #define DEBUG_SHUFFLE_BUFFER		1		// debug memory overruns with the shuffle buffer
 #define DEBUG_TICK_RANGES			0		// debug tick range parsing and evaluation
+#define DEBUG_TRAIT_DEMAND			0		// enable debugging logs about the trait "demand" evalutation process
 #define DEBUG_LESS_INTENSIVE		0		// decrease the frequency of some very intensive DEBUG checks
 
 
@@ -581,6 +582,10 @@ enum class TraitType : uint8_t {
 	kAdditive
 };
 
+std::string StringForTraitType(TraitType p_trait_type);
+TraitType TraitTypeForString(std::string type);				// raises if no match
+std::ostream& operator<<(std::ostream& p_out, TraitType p_trait_type);
+
 // This enumeration represents the type of a chromosome.  Note that the sex of an individual cannot always be inferred
 // from chromosomal state, and the user is allowed to play games with null haplosomes; the chromosomes follow the sex
 // of the individual, the sex of the individual does not follow the chromosomes.  See the initializeChromosome() doc.
@@ -1213,6 +1218,8 @@ extern const std::string &gStr_setSuppliedValue;
 extern const std::string &gStr_willAutolog;
 extern const std::string &gStr_context;
 
+extern const std::string gStr_additive;	// these trait type strings are not registered, no need
+extern const std::string gStr_multiplicative;
 extern const std::string gStr_A;	// these nucleotide strings are not registered, no need
 extern const std::string gStr_C;
 extern const std::string gStr_G;

core/slim_globals.h

@@ -1593,6 +1593,120 @@ late() { sim.killIndividuals(p1.subsetIndividuals(minAge=1)); }
 
 	SLiMAssertScriptSuccess(test_zygosity4);
 
+	
+	// Complex multi-trait, multi-chrom models that will (hopefully) exercise all the different code paths for
+	// phenotype evaluation -- callbacks, different DES combinations, neutral and non-neutral, etc.
+	// The intention here is that these models don't test themselves; they are stochastic and there is no
+	// expectation regarding their outcome.  Instead, _CheckPhenotypeForTrait() cross-checks them under DEBUG.
+	
+	std::string complex_multi_1 =	// this is an abbreviated version of test script complex_multi_test_1.slim
+	R"V0G0N(
+initialize() {
+	defineConstant("I1", 5.0);
+	defineConstant("I2", -5.0);
+	defineConstant("OPT1", 10.0);
+	defineConstant("OPT2", 10.0);
+	defineConstant("SD1", 2.0);
+	defineConstant("SD2", 2.0);
+	
+	initializeSex();
+	
+	// multiplicative traits
+	popgen1T = initializeTrait("popgen1T", "mul", 1.01, 1.0, 0.01, directFitnessEffect=T);  // will have a mix of dominance
+	popgen2T = initializeTrait("popgen2T", "mul", 1.01, 1.0, 0.01, directFitnessEffect=T);  // will be independent dominance
+	n1T = initializeTrait("n1T", "mul", NULL, NULL, NULL, directFitnessEffect=T);           // neutral with direct effect
+	n2T = initializeTrait("n2T", "mul", NULL, NULL, NULL, directFitnessEffect=F);           // neutral with no direct effect
+	
+	// additive traits
+	quant1T = initializeTrait("quant1T", "add", I1, 0.0, 0.01, directFitnessEffect=F);      // will have a mix of dominance
+	quant2T = initializeTrait("quant2T", "add", I2, 0.0, 0.01, directFitnessEffect=F);      // will be independent dominance
+	n3T = initializeTrait("n3T", "add", NULL, NULL, NULL, directFitnessEffect=F);           // non-neutral with no direct effect
+	
+	// quant1T / quant2T will be demanded in script; popgen1T / popgen2T / n1T will be demanded because they have direct effects
+	// calculation of popgen2T and quant2T should be extremely efficient since they are independent dominance
+	// calculation of n1T should be omitted entirely; SLiM should detect that it is neutral, and not even set phenotype values
+	// n2T and n3T should not be demanded, and should thus never be calculated by SLiM, which we can check in script
+	
+	// mutation types
+	initializeMutationType("m1", 0.4, "f", 0.0);                                  // neutral for all traits
+	
+	initializeMutationType("m2", 0.4, "e", 0.001);                                // beneficial for the popgen traits
+	m2.setEffectDistributionForTrait(c(n1T, n2T), "f", 0.0);                 		// neutral DES for the neutral traits
+	m2.setEffectDistributionForTrait(c(quant1T, quant2T), "n", 0.0, 0.1);         // unbiased normal DES for the additive traits
+	
+	initializeMutationType("m3", 0.4, "g", -0.001, 1.0);                          // deleterious for the popgen traits
+	m3.setEffectDistributionForTrait(c(n1T, n2T), "f", 0.0);                 		// neutral DES for the neutral traits
+	m3.setEffectDistributionForTrait(c(quant1T, quant2T), "n", 0.0, 0.1);         // unbiased normal DES for the additive traits
+	
+	c(m2,m3).setEffectDistributionForTrait(n3T, "n", -5.0, 0.5);                  // very biased and wide DES for n3T
+	
+	// set up independent dominance for popgen2T and quant2T; note that setting this for m1 should be unnecessary (it is neutral)
+	c(m1,m2,m3).setDefaultDominanceForTrait(c(popgen2T, quant2T), NAN);
+	
+	// prevent converting to substitutions for m2 and m3; once shifting the baseline offset is supported, this will not be needed
+	c(m2,m3).convertToSubstitution = F;
+	
+	initializeGenomicElementType("g1", m1, 1.0);          // neutral
+	initializeGenomicElementType("g2", 1:3, c(2, 1, 1));  // mixture
+	
+	ids = 1:5;
+	symbols = c(1, 2, "X", "Y", "MT");
+	lengths = rdunif(5, 1e7, 2e7);
+	types = c("A", "A", "X", "Y", "H");
+	names = c("A1", "A2", "X", "Y", "MT");
+	
+	for (id in ids, symbol in symbols, length in lengths, type in types, name in names)
+	{
+		initializeChromosome(id, length, type, symbol, name);
+		initializeMutationRate(1e-7);
+		initializeRecombinationRate(1e-8);
+		
+		if (id == 1)
+			initializeGenomicElement(g1);  // autosome 1 is pure neutral, using only m1
+		else
+			initializeGenomicElement(g2);  // autosome 2 is a mix, using m1 / m2 / m3
+	}
+}
+
+// set random dominance effects for the popgen1T and quant1T traits
+// other effects are generated as specified by the mutation type DES
+mutation(m2) {
+	mut.popgen1TDominance = runif(1);
+	mut.quant1TDominance = runif(1);
+	return T;
+}
+mutation(m3) {
+	mut.popgen1TDominance = runif(1);
+	mut.quant1TDominance = runif(1);
+	return T;
+}
+
+1 late() {
+	sim.addSubpop("p1", 20);
+}
+
+1: late() {
+	inds = sim.subpopulations.individuals;
+	sim.demandPhenotype(NULL, c(sim.quant1T, sim.quant2T));
+	fitnessEffect_q1 = dnorm(inds.quant1T, OPT1, SD1) / dnorm(0.0, 0.0, SD1);
+	fitnessEffect_q2 = dnorm(inds.quant2T, OPT2, SD2) / dnorm(0.0, 0.0, SD2);
+	inds.fitnessScaling = fitnessEffect_q1 * fitnessEffect_q2;
+}
+
+2: first() {
+	inds = sim.subpopulations.individuals;
+	
+	// check that traits that do not require calculation remain uncalculated
+	//if (!all(isNAN(inds.n1T))) stop("n1T was calculated unnecessarily");	// the smarts for this are not yet implemented!
+	if (!all(isNAN(inds.n2T))) stop("n2T was calculated unnecessarily");
+	if (!all(isNAN(inds.n3T))) stop("n3T was calculated unnecessarily");
+}
+
+50 late() { }
+	)V0G0N";
+	
+	SLiMAssertScriptSuccess(complex_multi_1);
+	
 	std::cout << "_RunMultitraitTests() done" << std::endl;
 }
 

core/slim_test_genetics.cpp

@@ -392,7 +392,7 @@ void Species::CheckOptimizationFlags(void)
 					slim_trait_index_t trait_index = trait->Index();
 					MutationTraitInfo &trait_info = mut_trait_info[trait_index];
 
-					if (trait_info.effect_size_ != 0.0)
+					if (trait_info.effect_size_ != (slim_effect_t)0.0)
 					{
 						// this mutation is non-neutral for this trait
 						if (trait->trait_all_neutral_mutations_ != false)
@@ -427,7 +427,7 @@ void Species::_NoteNonNeutralMutation(const Mutation *p_mut)
 		slim_trait_index_t trait_index = trait->Index();
 		MutationTraitInfo &trait_info = mut_trait_info[trait_index];
 
-		if (trait_info.effect_size_ != 0.0)
+		if (trait_info.effect_size_ != (slim_effect_t)0.0)
 		{
 			// this mutation is non-neutral for this trait
 			trait->trait_all_neutral_mutations_ = false;

core/species.cpp

@@ -1678,9 +1678,9 @@ EidosValue_SP Species::ExecuteContextFunction_initializeTrait(const std::string
 	std::string type_string = type_value->StringAtIndex_NOCAST(0, nullptr);
 	TraitType type;
 
-	if ((type_string == "multiplicative") || (type_string == "mul"))
+	if ((type_string == gStr_multiplicative) || (type_string == "mul"))
 		type = TraitType::kMultiplicative;
-	else if ((type_string == "additive") || (type_string == "add"))
+	else if ((type_string == gStr_additive) || (type_string == "add"))
 		type = TraitType::kAdditive;
 	else
 		EIDOS_TERMINATION << "ERROR (Species::ExecuteContextFunction_initializeTrait): initializeTrait() requires type to be either 'multiplicative' (or 'mul'), or 'additive' ('add')." << EidosTerminate();
@@ -3117,13 +3117,14 @@ EidosValue_SP Species::ExecuteMethod_demandPhenotype(EidosGlobalStringID p_metho
 	eidos_logical_t forceRecalc = forceRecalc_value->LogicalAtIndex_NOCAST(0, nullptr);
 
 	// subpops
-	std::vector<Subpopulation*> subpops_to_demand;
+	THREAD_SAFETY_IN_ACTIVE_PARALLEL("Species::ExecuteMethod_demandPhenotype(): usage of statics");
+	
+	static std::vector<Subpopulation*> subpops_to_demand;	// using and clearing a static prevents allocation thrash; should be safe from re-entry
+	subpops_to_demand.resize(0);
 
 	if (subpops_value->Type() == EidosValueType::kValueNULL)
 	{
 		// demand the specified phenotypes across all subpopulations
-		subpops_to_demand.resize(population_.subpops_.size());
-		
 		for (const std::pair<const slim_objectid_t,Subpopulation*> &subpop_pair : population_.subpops_)
 			subpops_to_demand.push_back(subpop_pair.second);
 	}
@@ -3134,17 +3135,25 @@ EidosValue_SP Species::ExecuteMethod_demandPhenotype(EidosGlobalStringID p_metho
 
 		if (requested_subpop_count)
 		{
-			subpops_to_demand.resize(requested_subpop_count);
-			
 			for (int requested_subpop_index = 0; requested_subpop_index < requested_subpop_count; ++requested_subpop_index)
-				subpops_to_demand.emplace_back(SLiM_ExtractSubpopulationFromEidosValue_io(subpops_value, requested_subpop_index, &community_, this, "demandPhenotype()"));		// SPECIES CONSISTENCY CHECK
+				subpops_to_demand.push_back(SLiM_ExtractSubpopulationFromEidosValue_io(subpops_value, requested_subpop_index, &community_, this, "demandPhenotype()"));		// SPECIES CONSISTENCY CHECK
 
 			// unique subpops_to_demand to avoid duplicated work
 			std::sort(subpops_to_demand.begin(), subpops_to_demand.end());
 			subpops_to_demand.resize(static_cast<size_t>(std::distance(subpops_to_demand.begin(), std::unique(subpops_to_demand.begin(), subpops_to_demand.end()))));
 		}
 	}
 
+#if DEBUG_TRAIT_DEMAND
+	std::cout << "# " << community_.Tick() << " --- demandPhenotype(): for traits {";
+	for (slim_trait_index_t trait_index : trait_indices)
+		std::cout << " " << Traits()[trait_index]->Name();
+	std::cout << " } in subpops {";
+	for (Subpopulation *subpop : subpops_to_demand)
+		std::cout << " p" << subpop->subpopulation_id_;
+	std::cout << " }, forceRecalc == " << (forceRecalc ? "T" : "F") << std::endl;
+#endif
+	
 	// validate non-neutral caches and independent-dominance precalculated values
 	// FIXME MULTITRAIT: VALIDATE NON-NEUTRAL CACHES HERE
 

core/species_eidos.cpp

@@ -1387,6 +1387,9 @@ void Subpopulation::UpdateFitness(std::vector<SLiMEidosBlock*> &p_subpop_mutatio
 	{
 		// we know this subpopulation has effectively constant fitness; we therefore don't express demand for
 		// any traits, which means trait may keep NAN values even if the traits have a direct fitness effect
+#if DEBUG_TRAIT_DEMAND
+		std::cout << "# " << community_.Tick() << " --- UpdateFitness() determined constant fitness of " << constant_fitness_value << std::endl;
+#endif
 
 		if (model_type_ == SLiMModelType::kModelTypeWF)
 		{
@@ -1425,11 +1428,24 @@ void Subpopulation::UpdateFitness(std::vector<SLiMEidosBlock*> &p_subpop_mutatio
 		// demand phenotypes for all the relevant traits
 		if (p_direct_effect_trait_indices.size())
 		{
+#if DEBUG_TRAIT_DEMAND
+				std::cout << "# " << community_.Tick() << " --- UpdateFitness() demanding traits {";
+				for (slim_trait_index_t trait_index : p_direct_effect_trait_indices)
+					std::cout << " " << species_.Traits()[trait_index]->Name();
+				std::cout << " } in subpop p" << subpopulation_id_  << ", forceRecalc == " << (p_force_trait_recalculation ? "T" : "F") << std::endl;
+#endif
+			
 			if (p_force_trait_recalculation)
 				Individual_Class::DemandPhenotype_SUBPOP<true>(&species_, this, p_direct_effect_trait_indices, p_subpop_mutationEffect_callbacks);
 			else
 				Individual_Class::DemandPhenotype_SUBPOP<false>(&species_, this, p_direct_effect_trait_indices, p_subpop_mutationEffect_callbacks);
 		}
+		else
+		{
+#if DEBUG_TRAIT_DEMAND
+				std::cout << "# " << community_.Tick() << " --- UpdateFitness() demanding NO traits in subpop p" << subpopulation_id_ << std::endl;
+#endif
+		}
 
 		// then loop over individuals and pull together the relevant phenotype values, fitnessEffect() callbacks,
 		// subpopulation fitnessScaling, and individual fitnessScaling to produce final individual fitness values;

core/subpopulation.cpp

@@ -113,8 +113,8 @@ EidosValue_SP Trait::GetProperty(EidosGlobalStringID p_property_id)
 			{
 				if (!static_type_string_multiplicative)
 				{
-					static_type_string_multiplicative = EidosValue_SP(new (gEidosValuePool->AllocateChunk()) EidosValue_String("multiplicative"));
-					static_type_string_additive = EidosValue_SP(new (gEidosValuePool->AllocateChunk()) EidosValue_String("additive"));
+					static_type_string_multiplicative = EidosValue_SP(new (gEidosValuePool->AllocateChunk()) EidosValue_String(gStr_multiplicative));
+					static_type_string_additive = EidosValue_SP(new (gEidosValuePool->AllocateChunk()) EidosValue_String(gStr_additive));
 				}
 			}