Multitrait metadata

VERSIONS

@@ -183,6 +183,13 @@ multitrait branch:
 		whenever a mutationEffect() callback is added, removed, rescheduled, or activated/deactivated in a way that affects the current tick (invalidating all affected trait values)
 		whenever the tick counter advances such that a previously active mutationEffect() callback becomes inactive, or vice versa (invalidating all affected trait values)
 		whenever the tick counter is changed in script to a new arbitrary value (invalidating all trait values in all individuals in all species)
+	tree-sequence recording work for multitrait:
+		individual metadata now has a variable-length per-trait array of phenotype and offset information
+		mutation metadata is now moved from the mutation table's metadata column to a table placed in top-level metadata with the json+struct codec
+		mutation metadata now has a variable-length per-trait array of effect size and dominance information
+		top-level JSON metadata now has a new `traits` key (required), which lists the traits used in the tree sequence; duplicated in provenance
+		the top-level JSON, mutation, and individual metadata schemas are updated accordingly
+		pushed the .trees file version from 0.9 to 1.0
 
 
 version 5.1 (Eidos version 4.1):

core/haplosome.cpp

@@ -4320,6 +4320,10 @@ EidosValue_SP Haplosome_Class::ExecuteMethod_removeMutations(EidosGlobalStringID
 					slim_position_t pos = mut->position_;
 
 					species->RecordNewDerivedState(target_haplosome, pos, empty_mut_vector);
+
+					// BCH 2/14/2026: All mutations removed need to be retained by the tree sequence until
+					// simplified away, so we can persist their metadata
+					species->NotifyMutationRemoved(mut);
 				}
 			}
 		}
@@ -4496,8 +4500,15 @@ EidosValue_SP Haplosome_Class::ExecuteMethod_removeMutations(EidosGlobalStringID
 				// TREE SEQUENCE RECORDING
 				// When doing tree recording, we additionally keep all fixed mutations (their ids) in a multimap indexed by their position
 				// This allows us to find all the fixed mutations at a given position quickly and easily, for calculating derived states
+				// BCH 2/14/2026: We now also need to remove the original mutation from our retained mutations list if it is there
 				if (species->RecordingTreeSequence())
+				{
 					pop.treeseq_substitutions_map_.emplace(mut->position_, sub);
+
+					// We just clear the flag here, and it is as if the mutation was actually removed from muts_retained_by_treeseq_;
+					// it will be garbage-collected the next time simplification occurs, so we don't have to search for it here
+					mut->retained_by_treeseq_ = false;
+				}
 
 				pop.substitutions_.emplace_back(sub);
 			}

core/individual.cpp

@@ -6395,10 +6395,10 @@ EidosValue_SP Individual_Class::ExecuteMethod_zygosityOfMutations(EidosGlobalStr
 	int64_t *integer_result_data = integer_result->data_mutable();
 
 	// tabulate the results for each individuals
-	const int8_t PRESENT_HETEROZYGOUS = 1;
-	const int8_t PRESENT_HOMOZYGOUS = 2;
-	const int8_t PRESENT_HEMIZYGOUS = 3;
-	const int8_t PRESENT_HAPLOID = 4;
+	const uint8_t PRESENT_HETEROZYGOUS = 1;
+	const uint8_t PRESENT_HOMOZYGOUS = 2;
+	const uint8_t PRESENT_HEMIZYGOUS = 3;
+	const uint8_t PRESENT_HAPLOID = 4;
 
 	for (int target_index = 0; target_index < target_size; ++target_index)
 	{

core/mutation.cpp

@@ -42,7 +42,7 @@ slim_mutationid_t gSLiM_next_mutation_id = 0;
 
 // This constructor is used when making a new mutation with effects and dominances provided by the caller; FIXME MULTITRAIT: needs to take a whole vector of each, per trait!
 Mutation::Mutation(MutationType *p_mutation_type_ptr, slim_chromosome_index_t p_chromosome_index, slim_position_t p_position, slim_effect_t p_selection_coeff, slim_effect_t p_dominance_coeff, slim_objectid_t p_subpop_index, slim_tick_t p_tick, int8_t p_nucleotide) :
-mutation_type_ptr_(p_mutation_type_ptr), position_(p_position), subpop_index_(p_subpop_index), origin_tick_(p_tick), chromosome_index_(p_chromosome_index), state_(MutationState::kNewMutation), nucleotide_(p_nucleotide), mutation_id_(gSLiM_next_mutation_id++)
+mutation_type_ptr_(p_mutation_type_ptr), position_(p_position), subpop_index_(p_subpop_index), origin_tick_(p_tick), chromosome_index_(p_chromosome_index), state_(MutationState::kNewMutation), retained_by_treeseq_(false), nucleotide_(p_nucleotide), mutation_id_(gSLiM_next_mutation_id++)
 {
 #ifdef DEBUG_LOCKS_ENABLED
 	mutation_block_LOCK.start_critical(2);
@@ -152,7 +152,7 @@ mutation_type_ptr_(p_mutation_type_ptr), position_(p_position), subpop_index_(p_
 
 // This constructor is used when making a new mutation with effects drawn from each trait's DES, and dominance taken from each trait's default dominance coefficient, both from the given mutation type
 Mutation::Mutation(MutationType *p_mutation_type_ptr, slim_chromosome_index_t p_chromosome_index, slim_position_t p_position, slim_objectid_t p_subpop_index, slim_tick_t p_tick, int8_t p_nucleotide) :
-mutation_type_ptr_(p_mutation_type_ptr), position_(p_position), subpop_index_(p_subpop_index), origin_tick_(p_tick), chromosome_index_(p_chromosome_index), state_(MutationState::kNewMutation), nucleotide_(p_nucleotide), mutation_id_(gSLiM_next_mutation_id++)
+mutation_type_ptr_(p_mutation_type_ptr), position_(p_position), subpop_index_(p_subpop_index), origin_tick_(p_tick), chromosome_index_(p_chromosome_index), state_(MutationState::kNewMutation), retained_by_treeseq_(false), nucleotide_(p_nucleotide), mutation_id_(gSLiM_next_mutation_id++)
 {
 #ifdef DEBUG_LOCKS_ENABLED
 	mutation_block_LOCK.start_critical(2);
@@ -295,9 +295,119 @@ mutation_type_ptr_(p_mutation_type_ptr), position_(p_position), subpop_index_(p_
 #endif
 }
 
+// This constructor is used when making a new mutation with effects and dominances provided by the caller, *and* a mutation id provided by the caller
+Mutation::Mutation(slim_mutationid_t p_mutation_id, MutationType *p_mutation_type_ptr, slim_chromosome_index_t p_chromosome_index, slim_position_t p_position, MutationTableMetadataRec *p_metadata_ptr) :
+mutation_type_ptr_(p_mutation_type_ptr), position_(p_position), subpop_index_(p_metadata_ptr->subpop_index_), origin_tick_(p_metadata_ptr->origin_tick_), chromosome_index_(p_chromosome_index), state_(MutationState::kNewMutation), retained_by_treeseq_(false), nucleotide_(p_metadata_ptr->nucleotide_), mutation_id_(p_mutation_id)
+{
+	Species &species = mutation_type_ptr_->species_;
+	const std::vector<Trait *> &traits = species.Traits();
+	MutationBlock *mutation_block = species.SpeciesMutationBlock();
+
+	// initialize the tag to the "unset" value
+	tag_value_ = SLIM_TAG_UNSET_VALUE;
+
+	// zero out our refcount and per-trait information, which is now kept in a separate buffer
+	MutationIndex mutation_index = mutation_block->IndexInBlock(this);
+	mutation_block->refcount_buffer_[mutation_index] = 0;
+
+	slim_trait_index_t trait_count = mutation_block->trait_count_;
+	MutationTraitInfo *mut_trait_info = mutation_block->TraitInfoForIndex(mutation_index);
+
+	// Below basically does the work of calling SetEffectSize() and SetDominance(), more efficiently since
+	// this is critical path.  See those methods for more comments on what is happening here.
+
+	is_neutral_for_all_traits_ = true;		// will be set to false below as needed
+	is_neutral_for_direct_fitness_traits_ = true;
+
+	// a dominance coefficient of NAN indicates independent dominance
+	independent_dominance_for_all_traits_ = true;
+	independent_dominance_for_any_traits_ = false;
+
+	for (slim_trait_index_t trait_index = 0; trait_index < trait_count; ++trait_index)
+	{
+		MutationTraitInfo *traitInfoRec = mut_trait_info + trait_index;
+		Trait *trait = traits[trait_index];
+		TraitType traitType = trait->Type();
+
+		// FIXME MULTITRAIT: This constructor needs to change to have a whole vector of trait information passed in, for effect and dominance
+		// for now we use the values passed in for trait 0, and make other traits neutral
+		slim_effect_t effect_size = p_metadata_ptr->per_trait_[trait_index].effect_size_;
+		slim_effect_t dominance = p_metadata_ptr->per_trait_[trait_index].dominance_coeff_;		// can be NAN
+		slim_effect_t hemizygous_dominance = p_metadata_ptr->per_trait_[trait_index].hemizygous_dominance_coeff_;
+
+		traitInfoRec->effect_size_ = effect_size;
+		traitInfoRec->dominance_coeff_UNSAFE_ = dominance;		// can be NAN
+		traitInfoRec->hemizygous_dominance_coeff_ = hemizygous_dominance;
+
+		if (effect_size != (slim_effect_t)0.0)
+		{
+			is_neutral_for_all_traits_ = false;
+
+			if (trait->HasDirectFitnessEffect())
+				is_neutral_for_direct_fitness_traits_ = false;
+
+			if (std::isnan(dominance))
+				independent_dominance_for_any_traits_ = true;
+			else
+				independent_dominance_for_all_traits_ = false;
+
+			// get the realized dominance to handle the possibility of independent dominance
+			slim_effect_t realized_dominance = RealizedDominanceForTrait(trait);
+
+			if (traitType == TraitType::kMultiplicative)
+			{
+				traitInfoRec->homozygous_effect_ = std::max((slim_effect_t)0.0, (slim_effect_t)1.0 + effect_size);
+				traitInfoRec->heterozygous_effect_ = std::max((slim_effect_t)0.0, (slim_effect_t)1.0 + realized_dominance * effect_size);
+				traitInfoRec->hemizygous_effect_ = std::max((slim_effect_t)0.0, (slim_effect_t)1.0 + hemizygous_dominance * effect_size);
+			}
+			else	// (traitType == TraitType::kAdditive)
+			{
+				traitInfoRec->homozygous_effect_ = ((slim_effect_t)2.0 * effect_size);
+				traitInfoRec->heterozygous_effect_ = ((slim_effect_t)2.0 * realized_dominance * effect_size);
+				traitInfoRec->hemizygous_effect_ = ((slim_effect_t)2.0 * hemizygous_dominance * effect_size);
+			}
+		}
+		else	// (effect == 0.0)
+		{
+			if (traitType == TraitType::kMultiplicative)
+			{
+				traitInfoRec->homozygous_effect_ = (slim_effect_t)1.0;
+				traitInfoRec->heterozygous_effect_ = (slim_effect_t)1.0;
+				traitInfoRec->hemizygous_effect_ = (slim_effect_t)1.0;
+			}
+			else	// (traitType == TraitType::kAdditive)
+			{
+				traitInfoRec->homozygous_effect_ = (slim_effect_t)0.0;
+				traitInfoRec->heterozygous_effect_ = (slim_effect_t)0.0;
+				traitInfoRec->hemizygous_effect_ = (slim_effect_t)0.0;
+			}
+		}
+	}
+
+	// this mutation will be added to the simulation somewhere, so tell the species about it
+	// (OK, it might not get added due to stacking policy or mutation() callbacks, but we assume it will be)
+	species.NoteChangedMutation(this);
+
+#if DEBUG
+	SelfConsistencyCheck(" in Mutation::Mutation()");
+#endif
+
+#if DEBUG_MUTATIONS()
+	std::cout << "Mutation constructed: " << this << std::endl;
+#endif
+
+	// Since a mutation id was supplied by the caller, we need to ensure that subsequent mutation ids generated do not collide
+	// This constructor (unlike the other Mutation() constructor above) is presently never called multithreaded,
+	// so we just enforce that here.  If that changes, it should start using the debug lock to detect races, as above.
+	THREAD_SAFETY_IN_ACTIVE_PARALLEL("Mutation::Mutation(): gSLiM_next_mutation_id change");
+
+	if (gSLiM_next_mutation_id <= mutation_id_)
+		gSLiM_next_mutation_id = mutation_id_ + 1;
+}
+
 // This constructor is used when making a new mutation with effects and dominances provided by the caller, *and* a mutation id provided by the caller
 Mutation::Mutation(slim_mutationid_t p_mutation_id, MutationType *p_mutation_type_ptr, slim_chromosome_index_t p_chromosome_index, slim_position_t p_position, slim_effect_t p_selection_coeff, slim_effect_t p_dominance_coeff, slim_objectid_t p_subpop_index, slim_tick_t p_tick, int8_t p_nucleotide) :
-mutation_type_ptr_(p_mutation_type_ptr), position_(p_position), subpop_index_(p_subpop_index), origin_tick_(p_tick), chromosome_index_(p_chromosome_index), state_(MutationState::kNewMutation), nucleotide_(p_nucleotide), mutation_id_(p_mutation_id)
+mutation_type_ptr_(p_mutation_type_ptr), position_(p_position), subpop_index_(p_subpop_index), origin_tick_(p_tick), chromosome_index_(p_chromosome_index), state_(MutationState::kNewMutation), retained_by_treeseq_(false), nucleotide_(p_nucleotide), mutation_id_(p_mutation_id)
 {
 	Species &species = mutation_type_ptr_->species_;
 	const std::vector<Trait *> &traits = species.Traits();

core/mutation.h

@@ -35,6 +35,7 @@
 
 class MutationType;
 class Trait;
+struct MutationTableMetadataRec;
 
 
 class Mutation_Class;
@@ -116,6 +117,18 @@ class Mutation : public EidosDictionaryRetained
 	const slim_tick_t origin_tick_;						// tick in which the mutation arose
 	slim_chromosome_index_t chromosome_index_;			// the (uint8_t) index of this mutation's chromosome
 	int state_ : 4;										// see MutationState above; 4 bits so we can represent -1
+	unsigned int scratch_ : 4;							// temporary scratch space for use by algorithms; regard as volatile outside your own code block
+
+	// This flag indicates that the mutation has been retained by the species because it belonged to a haplosome
+	// that was remembered by treeSeqRememberIndividuals(), or because it was removed in script or by stacking
+	// policy (and might therefore still be referenced by the tree sequence).  This flag is used by Species for
+	// bookkeeping.  It is used only when tree-sequence recording is enabled.
+	unsigned int retained_by_treeseq_ : 1;
+
+	// NOTE: there are three bits free here
+
+
+	// OPTIMIZATION FLAGS
 
 	// is_neutral_for_all_traits_ is true if all mutation effects are 0.0 (note a callback might override this).
 	// The state of is_neutral_for_all_traits_ is updated to reflect the state of the mutation when it changes.
@@ -147,8 +160,8 @@ class Mutation : public EidosDictionaryRetained
 	// if all effects for a mutation are neutral, this flag will be false.
 	unsigned int independent_dominance_for_any_traits_ : 1;
 
+
 	int8_t nucleotide_;									// the nucleotide being kept: A=0, C=1, G=2, T=3.  -1 is used to indicate non-nucleotide-based.
-	int8_t scratch_;									// temporary scratch space for use by algorithms; regard as volatile outside your own code block
 	const slim_mutationid_t mutation_id_;				// a unique id for each mutation, used to track mutations
 	slim_usertag_t tag_value_;							// a user-defined tag value
 
@@ -158,7 +171,7 @@ class Mutation : public EidosDictionaryRetained
 #endif
 
 #if DEBUG
-	mutable slim_refcount_t refcount_CHECK_;					// scratch space for checking of parallel refcounting
+	mutable slim_refcount_t refcount_CHECK_;				// scratch space for checking of parallel refcounting
 #endif
 
 	Mutation(const Mutation&) = delete;					// no copying
@@ -172,6 +185,9 @@ class Mutation : public EidosDictionaryRetained
 	// FIXME MULTITRAIT: needs to take a whole vector of each, per trait!
 	Mutation(MutationType *p_mutation_type_ptr, slim_chromosome_index_t p_chromosome_index, slim_position_t p_position, slim_effect_t p_selection_coeff, slim_effect_t p_dominance_coeff, slim_objectid_t p_subpop_index, slim_tick_t p_tick, int8_t p_nucleotide);
 
+	// This constructor is used by tree-sequence reading in Species::__CreateMutationsFromTabulation()
+	Mutation(slim_mutationid_t p_mutation_id, MutationType *p_mutation_type_ptr, slim_chromosome_index_t p_chromosome_index, slim_position_t p_position, MutationTableMetadataRec *p_metadata_ptr);
+
 	// This constructor is used when making a new mutation with effects and dominances PROVIDED by the caller, AND a mutation id provided by the caller
 	// FIXME MULTITRAIT: needs to take a whole vector of each, per trait!
 	Mutation(slim_mutationid_t p_mutation_id, MutationType *p_mutation_type_ptr, slim_chromosome_index_t p_chromosome_index, slim_position_t p_position, slim_effect_t p_selection_coeff, slim_effect_t p_dominance_coeff, slim_objectid_t p_subpop_index, slim_tick_t p_tick, int8_t p_nucleotide);

core/mutation_block.h

@@ -35,6 +35,7 @@
 
 #include "mutation.h"
 
+class Species;
 class Mutation;
 class MutationRun;
 

core/mutation_run.cpp

@@ -440,6 +440,8 @@ bool MutationRun::_EnforceStackPolicyForAddition(Mutation *p_mut_block_ptr, slim
 				if ((mut_position == p_position) && (mut->mutation_type_ptr_->stack_group_ == p_stack_group))
 				{
 					// The current scan position is a mutation that needs to be removed, so scan forward to skip copying it backward
+					// BCH 2/14/2026: this mutation would still be present in the tree sequence, so it needs to be retained by the species
+					mut->mutation_type_ptr_->species_.NotifyMutationRemoved(mut);
 					continue;
 				}
 				else

core/population.cpp

@@ -7677,6 +7677,7 @@ void Population::RemoveAllFixedMutations(void)
 		{
 			// When doing tree recording, we additionally keep all fixed mutations (their ids) in a multimap indexed by their position
 			// This allows us to find all the fixed mutations at a given position quickly and easily, for calculating derived states
+			// BCH 2/14/2026: We now also need to remove the original mutation from our retained mutations list if it is there
 			for (int i = 0; i < fixed_mutation_accumulator_size; i++)
 			{
 				Mutation *mut_to_remove = mut_block_ptr + fixed_mutation_accumulator[i];
@@ -7685,6 +7686,10 @@ void Population::RemoveAllFixedMutations(void)
 				species_.DoBaselineAccumulationForSubstitution(sub);
 				treeseq_substitutions_map_.emplace(mut_to_remove->position_, sub);
 				substitutions_.emplace_back(sub);
+
+				// We just clear the flag here, and it is as if the mutation was actually removed from muts_retained_by_treeseq_;
+				// it will be garbage-collected the next time simplification occurs, so we don't have to search for it here
+				mut_to_remove->retained_by_treeseq_ = false;
 			}
 		}
 		else