add genomic element lookup methods to Chromosome

Author: bhaller

QtSLiM/help/SLiMHelpClasses.html

@@ -85,8 +85,8 @@ geneConversionEnabled => (logical$)\
 
 \f4\fs20 \cf0 All of the 
 \f3\fs18 GenomicElement
-\f4\fs20  objects that comprise the chromosome.
-\f5 \
+\f4\fs20  objects that comprise the chromosome\cf2 , in sorted order (not necessarily in the order in which they were defined).
+\f5 \cf0 \
 \pard\pardeftab720\li720\fi-446\ri720\sb180\sa60\partightenfactor0
 
 \f3\fs18 \cf2 \expnd0\expndtw0\kerning0
@@ -532,7 +532,48 @@ It is generally recommended that the parent individual be supplied to this metho
 \
 \pard\pardeftab720\li720\fi-446\ri720\sb180\sa60\partightenfactor0
 
-\f3\fs18 \cf2 \'96\'a0(integer$)setAncestralNucleotides(is\'a0sequence)\
+\f3\fs18 \cf2 \kerning1\expnd0\expndtw0 \'96\'a0(object<GenomicElement>)genomicElementForPosition(integer\'a0positions)\
+\pard\pardeftab720\li547\ri720\sb60\sa60\partightenfactor0
+
+\f4\fs20 \cf2 Returns a vector of 
+\f3\fs18 GenomicElement
+\f4\fs20  objects corresponding to the given vector 
+\f3\fs18 positions
+\f4\fs20 , which contains base positions along the chromosome.  If every position lies within a defined genomic element, the returned vector will have the same length as 
+\f3\fs18 positions
+\f4\fs20 , and will correspond one-to-one with it.  However, if a position in 
+\f3\fs18 positions
+\f4\fs20  is not within a genomic element, no 
+\f3\fs18 GenomicElement
+\f4\fs20  object will be present for it in the returned vector, and so the returned vector will no longer have the same length as 
+\f3\fs18 positions
+\f4\fs20 , and will no longer correspond one-to-one with it.  The method 
+\f3\fs18 hasGenomicElementForPosition()
+\f4\fs20  can be used to detect this circumstance.\
+\pard\pardeftab720\li720\fi-446\ri720\sb180\sa60\partightenfactor0
+
+\f3\fs18 \cf2 \'96\'a0(logical)hasGenomicElementForPosition(integer\'a0positions)\
+\pard\pardeftab720\li547\ri720\sb60\sa60\partightenfactor0
+
+\f4\fs20 \cf2 Returns a 
+\f3\fs18 logical
+\f4\fs20  vector corresponding to the given vector 
+\f3\fs18 positions
+\f4\fs20 , which contains base positions along the chromosome.  The returned vector will have the same length as 
+\f3\fs18 positions
+\f4\fs20 , and will correspond one-to-one with it, containing 
+\f3\fs18 T
+\f4\fs20  if the corresponding position lies inside a genomic element, or 
+\f3\fs18 F
+\f4\fs20  if it does not.  The method 
+\f3\fs18 genomicElementForPosition()
+\f4\fs20  can be used to look up the 
+\f3\fs18 GenomicElement
+\f4\fs20  objects themselves.\
+\pard\pardeftab720\li720\fi-446\ri720\sb180\sa60\partightenfactor0
+
+\f3\fs18 \cf2 \expnd0\expndtw0\kerning0
+\'96\'a0(integer$)setAncestralNucleotides(is\'a0sequence)\
 \pard\pardeftab720\li547\ri720\sb60\sa60\partightenfactor0
 
 \f4\fs20 \cf2 This method, which may be called only in nucleotide-based models, replaces the ancestral nucleotide sequence for the model.  The 

SLiMgui/SLiMHelpClasses.rtf

@@ -31,6 +31,8 @@ development head (in the master branch):
 		for more details see https://doc.qt.io/qt-6/linux.html, https://doc.qt.io/qt-6/windows.html, https://doc.qt.io/qt-6/macos.html
 		building with CMake will now default to Qt6 if CMake can find it, and fall back to Qt5 otherwise; should be automatic
 	add a "Use OpenGL for fast display" checkbox in the Preferences panel, making it possible to disable OpenGL when it doesn't work well (#462)
+	add genomicElementForPosition() and hasGenomicElementForPosition() methods on Chromosome
+	policy change: genomicElements property is now in sorted order (not in order of declaration)
 
 
 version 4.2.2 (Eidos version 3.2.2):

VERSIONS

@@ -170,6 +170,12 @@ void Chromosome::InitializeDraws(void)
 	if (genomic_elements_.size() == 0)
 		EIDOS_TERMINATION << "ERROR (Chromosome::InitializeDraws): empty chromosome." << EidosTerminate();
 
+	// BCH 8/28/2024: we now sort the genomic elements here; we need them to be sorted later on anyway
+	// I avoided doing this before because it changed the behavior, but it is hard to imagine anybody caring
+	// Since elements must be non-overlapping, we only need to sort by start position
+	std::sort(genomic_elements_.begin(), genomic_elements_.end(),
+			  [](GenomicElement *e1, GenomicElement *e2) { return e1->start_position_ < e2->start_position_; });
+	
 	// determine which case we are working with: separate recombination maps for the sexes, or one map
 	auto rec_rates_H_size = recombination_rates_H_.size();
 	auto rec_rates_M_size = recombination_rates_M_.size();
@@ -614,12 +620,7 @@ void Chromosome::_InitializeOneMutationMap(gsl_ran_discrete_t *&p_lookup, std::v
 	// The class we use to represent these constant-rate subregions is GESubrange, declared in chromosome.h.
 	p_subranges.clear();
 
-	// We need to work with a *sorted* genomic elements vector here.  We sort it internally here, rather than in
-	// genomic_elements_, so that genomic_elements_ reflects exactly what they user gave us (mostly for backward
-	// compatibility).
-	std::vector<GenomicElement *> sorted_ge_vec = genomic_elements_;
-	
-	std::sort(sorted_ge_vec.begin(), sorted_ge_vec.end(), [](GenomicElement *ge1, GenomicElement *ge2) {return ge1->start_position_ < ge2->start_position_;});
+	// We need to work with a *sorted* genomic elements vector here.  BCH 8/28/2024: That is now guaranteed.
 
 	// This code deals in two different currencies: *requested* mutation rates and *adjusted* mutation rates.
 	// This stems from the fact that, beginning in SLiM 3.5, we unique mutation positions as we're generating
@@ -638,7 +639,7 @@ void Chromosome::_InitializeOneMutationMap(gsl_ran_discrete_t *&p_lookup, std::v
 	unsigned int mutrange_index = 0;
 	slim_position_t end_of_previous_mutrange = -1;
 
-	for (GenomicElement *ge_ptr : sorted_ge_vec)
+	for (GenomicElement *ge_ptr : genomic_elements_)
 	{
 		GenomicElement &ge = *ge_ptr;
 
@@ -1858,14 +1859,16 @@ EidosValue_SP Chromosome::ExecuteInstanceMethod(EidosGlobalStringID p_method_id,
 {
 	switch (p_method_id)
 	{
-		case gID_ancestralNucleotides:		return ExecuteMethod_ancestralNucleotides(p_method_id, p_arguments, p_interpreter);
-		case gID_setAncestralNucleotides:	return ExecuteMethod_setAncestralNucleotides(p_method_id, p_arguments, p_interpreter);
-		case gID_setGeneConversion:			return ExecuteMethod_setGeneConversion(p_method_id, p_arguments, p_interpreter);
-		case gID_setHotspotMap:				return ExecuteMethod_setHotspotMap(p_method_id, p_arguments, p_interpreter);
-		case gID_setMutationRate:			return ExecuteMethod_setMutationRate(p_method_id, p_arguments, p_interpreter);
-		case gID_setRecombinationRate:		return ExecuteMethod_setRecombinationRate(p_method_id, p_arguments, p_interpreter);
-		case gID_drawBreakpoints:			return ExecuteMethod_drawBreakpoints(p_method_id, p_arguments, p_interpreter);
-		default:							return super::ExecuteInstanceMethod(p_method_id, p_arguments, p_interpreter);
+		case gID_ancestralNucleotides:			return ExecuteMethod_ancestralNucleotides(p_method_id, p_arguments, p_interpreter);
+		case gID_genomicElementForPosition:		return ExecuteMethod_genomicElementForPosition(p_method_id, p_arguments, p_interpreter);
+		case gID_hasGenomicElementForPosition:	return ExecuteMethod_hasGenomicElementForPosition(p_method_id, p_arguments, p_interpreter);
+		case gID_setAncestralNucleotides:		return ExecuteMethod_setAncestralNucleotides(p_method_id, p_arguments, p_interpreter);
+		case gID_setGeneConversion:				return ExecuteMethod_setGeneConversion(p_method_id, p_arguments, p_interpreter);
+		case gID_setHotspotMap:					return ExecuteMethod_setHotspotMap(p_method_id, p_arguments, p_interpreter);
+		case gID_setMutationRate:				return ExecuteMethod_setMutationRate(p_method_id, p_arguments, p_interpreter);
+		case gID_setRecombinationRate:			return ExecuteMethod_setRecombinationRate(p_method_id, p_arguments, p_interpreter);
+		case gID_drawBreakpoints:				return ExecuteMethod_drawBreakpoints(p_method_id, p_arguments, p_interpreter);
+		default:								return super::ExecuteInstanceMethod(p_method_id, p_arguments, p_interpreter);
 	}
 }
 
@@ -2017,6 +2020,65 @@ EidosValue_SP Chromosome::ExecuteMethod_drawBreakpoints(EidosGlobalStringID p_me
 		return EidosValue_SP(new (gEidosValuePool->AllocateChunk()) EidosValue_Int(all_breakpoints));
 }
 
+GenomicElement *Chromosome::ElementForPosition(slim_position_t pos)
+{
+	auto element_iter = std::lower_bound(genomic_elements_.begin(), genomic_elements_.end(), pos,
+		[](const GenomicElement *e, slim_position_t position) { return e->end_position_ < position; });
+	
+	if (element_iter == genomic_elements_.end())
+		return nullptr;
+	
+	GenomicElement *element = *element_iter;
+	
+	if (element->start_position_ > pos)
+		return nullptr;
+	
+	return element;
+}
+
+//	*********************	(object<GenomicElement>)genomicElementForPosition(integer positions)
+//
+EidosValue_SP Chromosome::ExecuteMethod_genomicElementForPosition(EidosGlobalStringID p_method_id, const std::vector<EidosValue_SP> &p_arguments, EidosInterpreter &p_interpreter)
+{
+#pragma unused (p_method_id, p_arguments, p_interpreter)
+	EidosValue *positions_value = p_arguments[0].get();
+	int positions_count = positions_value->Count();
+	EidosValue_Object_SP obj_result_SP = EidosValue_Object_SP(new (gEidosValuePool->AllocateChunk()) EidosValue_Object(gSLiM_GenomicElement_Class));
+	EidosValue_Object *obj_result = obj_result_SP->reserve(positions_count);
+	const int64_t *positions_data = positions_value->IntData();
+	
+	for (int pos_index = 0; pos_index < positions_count; ++pos_index)
+	{
+		int64_t pos = positions_data[pos_index];
+		GenomicElement *element = ElementForPosition(pos);
+		
+		if (element)
+			obj_result->push_object_element_no_check_NORR(element);
+	}
+	
+	return obj_result_SP;
+}
+
+//	*********************	(logical)hasGenomicElementForPosition(integer positions)
+//
+EidosValue_SP Chromosome::ExecuteMethod_hasGenomicElementForPosition(EidosGlobalStringID p_method_id, const std::vector<EidosValue_SP> &p_arguments, EidosInterpreter &p_interpreter)
+{
+#pragma unused (p_method_id, p_arguments, p_interpreter)
+	EidosValue *positions_value = p_arguments[0].get();
+	int positions_count = positions_value->Count();
+	EidosValue_Logical_SP logical_result_SP = EidosValue_Logical_SP(new (gEidosValuePool->AllocateChunk()) EidosValue_Logical());
+	EidosValue_Logical *logical_result = logical_result_SP->reserve(positions_count);
+	const int64_t *positions_data = positions_value->IntData();
+	
+	for (int pos_index = 0; pos_index < positions_count; ++pos_index)
+	{
+		int64_t pos = positions_data[pos_index];
+		logical_result->push_logical_no_check(ElementForPosition(pos) ? true : false);
+	}
+	
+	return logical_result_SP;
+}
+
 //	*********************	(integer$)setAncestralNucleotides(is sequence)
 //
 EidosValue_SP Chromosome::ExecuteMethod_setAncestralNucleotides(EidosGlobalStringID p_method_id, const std::vector<EidosValue_SP> &p_arguments, EidosInterpreter &p_interpreter)
@@ -2577,6 +2639,8 @@ const std::vector<EidosMethodSignature_CSP> *Chromosome_Class::Methods(void) con
 
 		methods->emplace_back((EidosInstanceMethodSignature *)(new EidosInstanceMethodSignature(gStr_ancestralNucleotides, kEidosValueMaskInt | kEidosValueMaskString))->AddInt_OSN(gEidosStr_start, gStaticEidosValueNULL)->AddInt_OSN(gEidosStr_end, gStaticEidosValueNULL)->AddString_OS("format", EidosValue_String_SP(new (gEidosValuePool->AllocateChunk()) EidosValue_String("string"))));
 		methods->emplace_back((EidosInstanceMethodSignature *)(new EidosInstanceMethodSignature(gStr_drawBreakpoints, kEidosValueMaskInt))->AddObject_OSN("parent", gSLiM_Individual_Class, gStaticEidosValueNULL)->AddInt_OSN("n", gStaticEidosValueNULL));
+		methods->emplace_back((EidosInstanceMethodSignature *)(new EidosInstanceMethodSignature(gStr_genomicElementForPosition, kEidosValueMaskObject, gSLiM_GenomicElement_Class))->AddInt("positions"));
+		methods->emplace_back((EidosInstanceMethodSignature *)(new EidosInstanceMethodSignature(gStr_hasGenomicElementForPosition, kEidosValueMaskLogical))->AddInt("positions"));
 		methods->emplace_back((EidosInstanceMethodSignature *)(new EidosInstanceMethodSignature(gStr_setAncestralNucleotides, kEidosValueMaskInt | kEidosValueMaskSingleton))->AddIntString("sequence"));
 		methods->emplace_back((EidosInstanceMethodSignature *)(new EidosInstanceMethodSignature(gStr_setGeneConversion, kEidosValueMaskVOID))->AddNumeric_S("nonCrossoverFraction")->AddNumeric_S("meanLength")->AddNumeric_S("simpleConversionFraction")->AddNumeric_OS("bias", gStaticEidosValue_Integer0));
 		methods->emplace_back((EidosInstanceMethodSignature *)(new EidosInstanceMethodSignature(gStr_setHotspotMap, kEidosValueMaskVOID))->AddNumeric("multipliers")->AddInt_ON("ends", gStaticEidosValueNULL)->AddString_OS("sex", gStaticEidosValue_StringAsterisk));

core/chromosome.cpp

@@ -60,6 +60,7 @@ class Chromosome : public EidosDictionaryRetained
 private:
 #endif
 
+	// This vector contains all the genomic elements for this chromosome.  It is in sorted order once initialization is complete.
 	std::vector<GenomicElement *> genomic_elements_;		// OWNED POINTERS: genomic elements belong to the chromosome
 
 	// We now allow different recombination maps for males and females, optionally.  Unfortunately, this means we have a bit of an
@@ -229,6 +230,9 @@ class Chromosome : public EidosDictionaryRetained
 												   double &p_both_0, double &p_both_0_OR_mut_0_break_non0, double &p_both_0_OR_mut_0_break_non0_OR_mut_non0_break_0);
 #endif
 
+	// looking up genomic elements quickly, by binary search
+	GenomicElement *ElementForPosition(slim_position_t pos);
+	
 	// internal methods for throwing errors from inline functions when assumptions about the configuration of maps are violated
 	void MutationMapConfigError(void) const __attribute__((__noreturn__)) __attribute__((cold)) __attribute__((analyzer_noreturn));
 	void RecombinationMapConfigError(void) const __attribute__((__noreturn__)) __attribute__((cold)) __attribute__((analyzer_noreturn));
@@ -250,6 +254,8 @@ class Chromosome : public EidosDictionaryRetained
 
 	virtual EidosValue_SP ExecuteInstanceMethod(EidosGlobalStringID p_method_id, const std::vector<EidosValue_SP> &p_arguments, EidosInterpreter &p_interpreter) override;
 	EidosValue_SP ExecuteMethod_ancestralNucleotides(EidosGlobalStringID p_method_id, const std::vector<EidosValue_SP> &p_arguments, EidosInterpreter &p_interpreter);
+	EidosValue_SP ExecuteMethod_genomicElementForPosition(EidosGlobalStringID p_method_id, const std::vector<EidosValue_SP> &p_arguments, EidosInterpreter &p_interpreter);
+	EidosValue_SP ExecuteMethod_hasGenomicElementForPosition(EidosGlobalStringID p_method_id, const std::vector<EidosValue_SP> &p_arguments, EidosInterpreter &p_interpreter);
 	EidosValue_SP ExecuteMethod_setAncestralNucleotides(EidosGlobalStringID p_method_id, const std::vector<EidosValue_SP> &p_arguments, EidosInterpreter &p_interpreter);
 	EidosValue_SP ExecuteMethod_setGeneConversion(EidosGlobalStringID p_method_id, const std::vector<EidosValue_SP> &p_arguments, EidosInterpreter &p_interpreter);
 	EidosValue_SP ExecuteMethod_setHotspotMap(EidosGlobalStringID p_method_id, const std::vector<EidosValue_SP> &p_arguments, EidosInterpreter &p_interpreter);

core/chromosome.h

@@ -1271,6 +1271,8 @@ const std::string &gStr_maxDistance = EidosRegisteredString("maxDistance", gID_m
 // mostly method names
 const std::string &gStr_ancestralNucleotides = EidosRegisteredString("ancestralNucleotides", gID_ancestralNucleotides);
 const std::string &gStr_nucleotides = EidosRegisteredString("nucleotides", gID_nucleotides);
+const std::string &gStr_genomicElementForPosition = EidosRegisteredString("genomicElementForPosition", gID_genomicElementForPosition);
+const std::string &gStr_hasGenomicElementForPosition = EidosRegisteredString("hasGenomicElementForPosition", gID_hasGenomicElementForPosition);
 const std::string &gStr_setAncestralNucleotides = EidosRegisteredString("setAncestralNucleotides", gID_setAncestralNucleotides);
 const std::string &gStr_setGeneConversion = EidosRegisteredString("setGeneConversion", gID_setGeneConversion);
 const std::string &gStr_setHotspotMap = EidosRegisteredString("setHotspotMap", gID_setHotspotMap);

core/slim_globals.cpp

@@ -865,6 +865,8 @@ extern const std::string &gStr_maxDistance;
 
 extern const std::string &gStr_ancestralNucleotides;
 extern const std::string &gStr_nucleotides;
+extern const std::string &gStr_genomicElementForPosition;
+extern const std::string &gStr_hasGenomicElementForPosition;
 extern const std::string &gStr_setAncestralNucleotides;
 extern const std::string &gStr_setGeneConversion;
 extern const std::string &gStr_setHotspotMap;
@@ -1283,6 +1285,8 @@ enum _SLiMGlobalStringID : int {
 
 	gID_ancestralNucleotides,
 	gID_nucleotides,
+	gID_genomicElementForPosition,
+	gID_hasGenomicElementForPosition,
 	gID_setAncestralNucleotides,
 	gID_setGeneConversion,
 	gID_setHotspotMap,