extend add[Multi]Recombinant() with breakpoint drawing, gene conversi…

…on support

VERSIONS

@@ -93,6 +93,9 @@ development head (in the master branch):
 		policy change: randomizeStrands must now usually be explicitly specified for both addRecombinant() and addMultiRecombinant();
 			the default is now NULL, and NULL errors so that T or F must be explicitly given unless it does not matter (i.e., there is no recombination)
 		implement addMultiRecombinant() for complete control over reproduction in multi-chromosome models
+		extend addRecombinant() and addMultiRecombinant() to allow NULL as a breaks vector with two strands, meaning "draw breakpoints for me", to avoid needing drawBreakpoints() in most cases
+			this new code path also allows gene conversion including heteroduplex mismatch repair and gBGC to work with addRecombinant() and addMultiRecombinant()
+			also change to allow a recombination rate to be set for haploid chromosomes; this is potentially useful with addRecombinant() etc.
 
 
 version 4.3 (Eidos version 3.3):

core/chromosome.h

@@ -306,7 +306,6 @@ class Chromosome : public EidosDictionaryRetained
 	inline bool UsingSingleMutationMap(void) const { return single_mutation_map_; }
 	inline size_t GenomicElementCount(void) const { return genomic_elements_.size(); }
 
-	bool RequiresZeroRecombination(void) const;
 	bool DefaultsToZeroRecombination(void) const;
 
 	// draw the number of mutations that occur, based on the overall mutation rate
@@ -325,9 +324,12 @@ class Chromosome : public EidosDictionaryRetained
 	// draw the number of breakpoints that occur, based on the overall recombination rate
 	int DrawBreakpointCount(IndividualSex p_sex) const;
 
-	// choose a set of recombination breakpoints, based on recomb. intervals, overall recomb. rate, and gene conversion parameters
-	void DrawCrossoverBreakpoints(IndividualSex p_parent_sex, const int p_num_breakpoints, std::vector<slim_position_t> &p_crossovers) const;
-	void DrawDSBBreakpoints(IndividualSex p_parent_sex, const int p_num_breakpoints, std::vector<slim_position_t> &p_crossovers, std::vector<slim_position_t> &p_heteroduplex) const;
+	// choose a set of recombination breakpoints, based on recomb. intervals, overall recomb. rate, and gene
+	// conversion parameters; DrawBreakpoints() is the high-level funnel that most callers should use, whereas
+	// the low-level functions do not handle recombination() callbacks and other niceties
+	void _DrawCrossoverBreakpoints(IndividualSex p_parent_sex, const int p_num_breakpoints, std::vector<slim_position_t> &p_crossovers) const;
+	void _DrawDSBBreakpoints(IndividualSex p_parent_sex, const int p_num_breakpoints, std::vector<slim_position_t> &p_crossovers, std::vector<slim_position_t> &p_heteroduplex) const;
+	void DrawBreakpoints(Individual *p_parent, Haplosome *p_haplosome1, Haplosome *p_haplosome2, int p_num_breakpoints, std::vector<slim_position_t> &p_crossovers, std::vector<slim_position_t> *p_heteroduplex, const char *p_caller_name);
 
 #ifndef USE_GSL_POISSON
 	// draw both the mutation count and breakpoint count, using a single Poisson draw for speed

core/population.cpp

@@ -2459,7 +2459,7 @@ void Population::EvolveSubpopulation(Subpopulation &p_subpop, bool p_mate_choice
 }
 
 // apply recombination() callbacks to a generated child; a return of true means breakpoints were changed
-bool Population::ApplyRecombinationCallbacks(Haplosome *p_haplosome1, Haplosome *p_haplosome2, std::vector<slim_position_t> &p_crossovers, std::vector<SLiMEidosBlock*> &p_recombination_callbacks)
+bool Population::ApplyRecombinationCallbacks(Individual *p_parent, Haplosome *p_haplosome1, Haplosome *p_haplosome2, std::vector<slim_position_t> &p_crossovers, std::vector<SLiMEidosBlock*> &p_recombination_callbacks)
 {
 	THREAD_SAFETY_IN_ANY_PARALLEL("Population::ApplyRecombinationCallbacks(): running Eidos callback");
 
@@ -2518,16 +2518,13 @@ bool Population::ApplyRecombinationCallbacks(Haplosome *p_haplosome1, Haplosome
 			// the value objects, and we know that the values we are setting here will not change (the objects
 			// referred to by the values may change, but the values themselves will not change).
 			if (recombination_callback->contains_individual_)
-			{
-				Individual *individual = p_haplosome1->OwningIndividual();
-				callback_symbols.InitializeConstantSymbolEntry(gID_individual, individual->CachedEidosValue());
-			}
+				callback_symbols.InitializeConstantSymbolEntry(gID_individual, p_parent->CachedEidosValue());
 			if (recombination_callback->contains_haplosome1_)
 				callback_symbols.InitializeConstantSymbolEntry(gID_haplosome1, p_haplosome1->CachedEidosValue());
 			if (recombination_callback->contains_haplosome2_)
 				callback_symbols.InitializeConstantSymbolEntry(gID_haplosome2, p_haplosome2->CachedEidosValue());
 			if (recombination_callback->contains_subpop_)
-				callback_symbols.InitializeConstantSymbolEntry(gID_subpop, p_haplosome1->OwningIndividual()->subpopulation_->SymbolTableEntry().second);
+				callback_symbols.InitializeConstantSymbolEntry(gID_subpop, p_parent->subpopulation_->SymbolTableEntry().second);
 
 			// All the variable entries for the crossovers and gene conversion start/end points
 			if (recombination_callback->contains_breakpoints_)
@@ -2649,8 +2646,9 @@ void Population::HaplosomeCrossed(Chromosome &p_chromosome, Haplosome &p_child_h
 
 	// some behaviors -- which callbacks to use, which recombination/mutation rate to use, subpop of origin for mutations, etc. --
 	// depend upon characteristics of the first parent, so we fetch the necessary properties here
-	Subpopulation *source_subpop = parent_haplosome_1->individual_->subpopulation_;
-	IndividualSex parent_sex = parent_haplosome_1->individual_->sex_;
+	Individual *parent_individual = parent_haplosome_1->individual_;
+	Subpopulation *source_subpop = parent_individual->subpopulation_;
+	IndividualSex parent_sex = parent_individual->sex_;
 
 	// determine how many mutations and breakpoints we have
 	int num_mutations, num_breakpoints;
@@ -2682,25 +2680,29 @@ void Population::HaplosomeCrossed(Chromosome &p_chromosome, Haplosome &p_child_h
 		//std::cout << num_mutations << " mutations, " << num_breakpoints << " breakpoints" << std::endl;
 
 		// draw the breakpoints based on the recombination rate map, and sort and unique the result
+		// we don't use Chromosome::DrawBreakpoints(), for speed, but this code mirrors it
 		if (num_breakpoints)
 		{
 			if (p_chromosome.using_DSB_model_)
-				p_chromosome.DrawDSBBreakpoints(parent_sex, num_breakpoints, all_breakpoints, heteroduplex);
+				p_chromosome._DrawDSBBreakpoints(parent_sex, num_breakpoints, all_breakpoints, heteroduplex);
 			else
-				p_chromosome.DrawCrossoverBreakpoints(parent_sex, num_breakpoints, all_breakpoints);
+				p_chromosome._DrawCrossoverBreakpoints(parent_sex, num_breakpoints, all_breakpoints);
+
+			// all_breakpoints is sorted and uniqued at this point
 
 			if (f_callbacks && p_recombination_callbacks)
 			{
 				// a non-zero number of breakpoints, with recombination callbacks
 				if (p_chromosome.using_DSB_model_ && (p_chromosome.simple_conversion_fraction_ != 1.0))
 					EIDOS_TERMINATION << "ERROR (Population::HaplosomeCrossed): recombination() callbacks may not be used when complex gene conversion tracts are in use, since recombination() callbacks have no support for heteroduplex regions." << EidosTerminate();
 
-				ApplyRecombinationCallbacks(parent_haplosome_1, parent_haplosome_2, all_breakpoints, *p_recombination_callbacks);
+				bool breaks_changed = ApplyRecombinationCallbacks(parent_individual, parent_haplosome_1, parent_haplosome_2, all_breakpoints, *p_recombination_callbacks);
 				num_breakpoints = (int)all_breakpoints.size();
 
 				if (num_breakpoints)
 				{
-					if (all_breakpoints.size() > 1)
+					// we only sort/unique if the breakpoints have changed, since they were sorted/uniqued before
+					if (breaks_changed && (all_breakpoints.size() > 1))
 					{
 						std::sort(all_breakpoints.begin(), all_breakpoints.end());
 						all_breakpoints.erase(unique(all_breakpoints.begin(), all_breakpoints.end()), all_breakpoints.end());
@@ -2728,7 +2730,7 @@ void Population::HaplosomeCrossed(Chromosome &p_chromosome, Haplosome &p_child_h
 			if (p_chromosome.using_DSB_model_ && (p_chromosome.simple_conversion_fraction_ != 1.0))
 				EIDOS_TERMINATION << "ERROR (Population::HaplosomeCrossed): recombination() callbacks may not be used when complex gene conversion tracts are in use, since recombination() callbacks have no support for heteroduplex regions." << EidosTerminate();
 
-			ApplyRecombinationCallbacks(parent_haplosome_1, parent_haplosome_2, all_breakpoints, *p_recombination_callbacks);
+			ApplyRecombinationCallbacks(parent_individual, parent_haplosome_1, parent_haplosome_2, all_breakpoints, *p_recombination_callbacks);
 			num_breakpoints = (int)all_breakpoints.size();
 
 			if (num_breakpoints)
@@ -2757,6 +2759,15 @@ void Population::HaplosomeCrossed(Chromosome &p_chromosome, Haplosome &p_child_h
 		}
 	}
 
+	// A leading zero in the breakpoints vector switches copy strands before copying begins.
+	// In HaplosomeRecombined() we explicitly remove that, to prevent it from being recorded
+	// in the tree sequence.  That's a bit trickier to do here, since we've deferred adding the
+	// past-the-end marker above in some cases.  It is not essential to do, I think, just a
+	// nicety that makes for a cleaner recorded tree sequence.  Also, a leading zero should
+	// be much less common in HaplosomeCrossed(), since the user is not playing games the way
+	// they do with HaplosomeRecombined().  So I'm going to skip it for now, which I think is
+	// basically harmless; but I'll revisit this later.  FIXME MULTICHROM
+
 	// TREE SEQUENCE RECORDING
 	bool recording_tree_sequence = f_treeseq;
 	bool recording_tree_sequence_mutations = f_treeseq && species_.RecordingTreeSequenceMutations();
@@ -3803,8 +3814,24 @@ void Population::HaplosomeRecombined(Chromosome &p_chromosome, Haplosome &p_chil
 	// determine how many mutations we have
 	int num_mutations = p_chromosome.DrawMutationCount(parent_sex);
 
-	// we need a defined end breakpoint, so we add it now
-	p_breakpoints.emplace_back(p_chromosome.last_position_mutrun_ + 1);
+	// A leading zero in the breakpoints vector switches copy strands before copying begins.
+	// We want to handle that up front, primarily because we don't want to record it in treeseq.
+	// FIXME MULTICHROM erase() here can be very expensive; would be good to switch to using a start pointer
+	// and count, which would allow the first element to be removed in O(1) time by advancing the pointer
+	// I've put off doing that change because it echos across many calls and will make big diffs
+	// As a part of this work, RecordNewHaplosome() should get a derived version for the no-recombination case,
+	// especially the null haplosome case; recording those should be fast to avoid bogging down the model
+	while (p_breakpoints.front() == 0)
+	{
+		p_breakpoints.erase(p_breakpoints.begin());
+		std::swap(parent_haplosome_1, parent_haplosome_2);
+	}
+
+	// we need a defined end breakpoint, so we add it now; we don't want more than one, though
+	// in some cases the caller might already have this breakpoint added, so we need to check
+	// (that happens with multiple children in addRecombinant(), for example)
+	if (p_breakpoints.back() <= p_chromosome.last_position_mutrun_)
+		p_breakpoints.emplace_back(p_chromosome.last_position_mutrun_ + 1);
 
 	// TREE SEQUENCE RECORDING
 	bool recording_tree_sequence = f_treeseq;

core/population.h

@@ -220,7 +220,8 @@ class Population
 	bool ApplyModifyChildCallbacks(Individual *p_child, Individual *p_parent1, Individual *p_parent2, bool p_is_selfing, bool p_is_cloning, Subpopulation *p_target_subpop, Subpopulation *p_source_subpop, std::vector<SLiMEidosBlock*> &p_modify_child_callbacks);
 
 	// apply recombination() callbacks to a generated child; a return of true means the breakpoints were changed
-	bool ApplyRecombinationCallbacks(Haplosome *p_haplosome1, Haplosome *p_haplosome2, std::vector<slim_position_t> &p_crossovers, std::vector<SLiMEidosBlock*> &p_recombination_callbacks);
+	// note that in the general case (e.g., addRecombinant()), the haplosomes might not belong to the parent
+	bool ApplyRecombinationCallbacks(Individual *p_parent, Haplosome *p_haplosome1, Haplosome *p_haplosome2, std::vector<slim_position_t> &p_crossovers, std::vector<SLiMEidosBlock*> &p_recombination_callbacks);
 
 	// generate a child haplosome from parental haplosome(s), very directly -- no null haplosomes etc., just cross/clone/recombine
 	// these methods are templated with variants for speed; see also MungeIndividualCrossed() etc.

core/species.cpp

@@ -5111,10 +5111,20 @@ void Species::RecordNewHaplosome(std::vector<slim_position_t> *p_breakpoints, Ha
 		right = (*p_breakpoints)[i];
 
 		tsk_id_t parent = (tsk_id_t) (polarity ? haplosome1TSKID : haplosome2TSKID);
+		polarity = !polarity;
+
+		// Sometimes the user might add a breakpoint at 0, to flip the initial copy strand, as in the meiotic
+		// drive recipe.  If they do that, a left==right breakpoint might make it in to here.  That would be
+		// a bug in the caller.  This has never been seen in the wild, so I'll make it DEBUG only.  In non-
+		// DEBUG runs the tree sequence will fail to pass integrity checks, with TSK_ERR_BAD_EDGE_INTERVAL.
+#if DEBUG
+		if (left >= right)
+			EIDOS_TERMINATION << "ERROR (Species::RecordNewHaplosome): (internal error) a left==right breakpoint was passed to RecordNewHaplosome()." << EidosTerminate();
+#endif
+
 		int ret = tsk_edge_table_add_row(&tsinfo.tables_.edges, left, right, parent, offspringTSKID, NULL, 0);
 		if (ret < 0) handle_error("tsk_edge_table_add_row", ret);
 
-		polarity = !polarity;
 		left = right;
 	}
 

core/species_eidos.cpp

@@ -667,13 +667,11 @@ EidosValue_SP Species::ExecuteContextFunction_initializeRecombinationRate(const
 
 		double recombination_rate = rates_value->NumericAtIndex_NOCAST(0, nullptr);
 
-		// check values
+		// check values; I thought about requiring a rate of 0.0 for all haploid chromosome types, but maybe
+		// the user wants to recombine them sometimes with addRecombinant(), no need to prevent them
 		if ((recombination_rate < 0.0) || (recombination_rate > 0.5) || std::isnan(recombination_rate))
 			EIDOS_TERMINATION << "ERROR (Species::ExecuteContextFunction_initializeRecombinationRate): initializeRecombinationRate() requires rates to be in [0.0, 0.5] (" << EidosStringForFloat(recombination_rate) << " supplied)." << EidosTerminate();
 
-		if ((recombination_rate != 0.0) && chromosome->RequiresZeroRecombination())
-			EIDOS_TERMINATION << "ERROR (Chromosome::ExecuteMethod_setRecombinationRate): setRecombinationRate() requires a recombination rate of 0.0 for all chromosome types except 'A', 'H', 'X', and 'Z'." << EidosTerminate();
-
 		// then adopt them
 		rates.clear();
 		positions.clear();
@@ -701,9 +699,6 @@ EidosValue_SP Species::ExecuteContextFunction_initializeRecombinationRate(const
 			if ((recombination_rate < 0.0) || (recombination_rate > 0.5) || std::isnan(recombination_rate))
 				EIDOS_TERMINATION << "ERROR (Species::ExecuteContextFunction_initializeRecombinationRate): initializeRecombinationRate() requires rates to be in [0.0, 0.5] (" << EidosStringForFloat(recombination_rate) << " supplied)." << EidosTerminate();
 
-			if ((recombination_rate != 0.0) && chromosome->RequiresZeroRecombination())
-				EIDOS_TERMINATION << "ERROR (Chromosome::ExecuteMethod_setRecombinationRate): setRecombinationRate() requires a recombination rate of 0.0 for all chromosome types except 'A', 'H', 'X', and 'Z'." << EidosTerminate();
-
 			if (chromosome->extent_immutable_)
 			{
 				if ((recombination_end_position <= 0) || (recombination_end_position > chromosome->last_position_))