Skip to content

Commit eef9e92

Browse files
committed
moved duplicate haplosome check to only alive individuals; closes #538
1 parent e0baffb commit eef9e92

1 file changed

Lines changed: 12 additions & 21 deletions

File tree

core/species.cpp

Lines changed: 12 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -8725,6 +8725,7 @@ void Species::__CreateSubpopulationsFromTabulation(std::unordered_map<slim_objec
87258725
// This could be done with a hash table, but I imagine that would be slower until the number of individuals becomes very large
87268726
// Also, I'm a bit nervous about putting a large number of consecutive integers into a hash table, re: edge-case performance
87278727
std::vector<slim_pedigreeid_t> pedigree_id_check;
8728+
std::vector<slim_haplosomeid_t> haplosome_id_check;
87288729

87298730
gSLiM_next_pedigree_id = 0;
87308731

@@ -8804,12 +8805,16 @@ void Species::__CreateSubpopulationsFromTabulation(std::unordered_map<slim_objec
88048805
individual->spatial_z_ = subpop_info.spatial_z_[tabulation_index];
88058806

88068807
p_nodeToHaplosomeMap.emplace(node_id_0, individual->haplosomes_[first_haplosome_index]);
8807-
individual->haplosomes_[first_haplosome_index]->haplosome_id_ = pedigree_id * 2;
8808+
slim_haplosomeid_t haplosome_id = pedigree_id * 2;
8809+
individual->haplosomes_[first_haplosome_index]->haplosome_id_ = haplosome_id;
8810+
haplosome_id_check.emplace_back(haplosome_id); // we will test for collisions below
88088811

88098812
if (last_haplosome_index != first_haplosome_index)
88108813
{
88118814
p_nodeToHaplosomeMap.emplace(node_id_1, individual->haplosomes_[last_haplosome_index]);
8812-
individual->haplosomes_[last_haplosome_index]->haplosome_id_ = pedigree_id * 2 + 1;
8815+
haplosome_id = pedigree_id * 2 + 1;
8816+
individual->haplosomes_[last_haplosome_index]->haplosome_id_ = haplosome_id;
8817+
haplosome_id_check.emplace_back(haplosome_id); // we will test for collisions below
88138818
}
88148819

88158820
// check the referenced nodes; right now this is not essential for re-creating the saved state, but is just a crosscheck
@@ -9812,13 +9817,12 @@ void Species::__CheckNodePedigreeIDs(__attribute__((unused)) EidosInterpreter *p
98129817
// could be nodes in the node table with haplosome pedigree IDs greater than those in use by individuals, in nonWF models.
98139818
// See https://github.com/MesserLab/SLiM/pull/420 for an example model that does this very easily.
98149819

9815-
// Also, check for duplicate pedigree IDs, just in case. __CreateSubpopulationsFromTabulation() does this for individual
9816-
// pedigree IDs; we do it for node pedigree IDs. I decided to use a vector with std::sort() to check even though it is
9817-
// O(n log n), rather than a hash table for O(n), because I'm nervous about hitting a bad edge case with the hash table
9818-
// due to the nature of the values being inserted. Shouldn't be a big deal in the grand scheme of things.
9820+
// Previously, we checked for duplicate pedigree IDs here as well, just in case.
9821+
// __CreateSubpopulationsFromTabulation() does this for haplosome IDs in
9822+
// living individuals already pedigree IDs; however, it was found to be
9823+
// overly restrictive, in situations involving merging of parallel simulations.
98199824
tsk_node_table_t &node_table = tables.nodes;
98209825
tsk_size_t node_count = node_table.num_rows;
9821-
std::vector<slim_haplosomeid_t> haplosome_id_check;
98229826

98239827
for (tsk_size_t j = 0; (size_t)j < node_count; j++)
98249828
{
@@ -9831,11 +9835,7 @@ void Species::__CheckNodePedigreeIDs(__attribute__((unused)) EidosInterpreter *p
98319835
{
98329836
// get the metadata record and check the haplosome pedigree ID
98339837
HaplosomeMetadataRec *metadata_rec = (HaplosomeMetadataRec *)(node_table.metadata + offset1);
9834-
slim_haplosomeid_t haplosome_id = metadata_rec->haplosome_id_;
9835-
9836-
haplosome_id_check.emplace_back(haplosome_id); // we will test for collisions below
9837-
9838-
slim_pedigreeid_t pedigree_id = haplosome_id / 2; // rounds down to integer
9838+
slim_pedigreeid_t pedigree_id = metadata_rec->haplosome_id_ / 2; // rounds down to integer
98399839

98409840
if (pedigree_id >= gSLiM_next_pedigree_id)
98419841
{
@@ -9852,15 +9852,6 @@ void Species::__CheckNodePedigreeIDs(__attribute__((unused)) EidosInterpreter *p
98529852
}
98539853
}
98549854
}
9855-
9856-
/* REMOVE TO SEE IF IT FIXES THE PROBLEM
9857-
// Check for haplosome pedigree ID collisions by sorting and looking for duplicates
9858-
std::sort(haplosome_id_check.begin(), haplosome_id_check.end());
9859-
const auto duplicate = std::adjacent_find(haplosome_id_check.begin(), haplosome_id_check.end());
9860-
9861-
if (duplicate != haplosome_id_check.end())
9862-
EIDOS_TERMINATION << "ERROR (Species::__CheckNodePedigreeIDs): the haplosome pedigree ID value " << *duplicate << " was used more than once; haplosome pedigree IDs must be unique." << EidosTerminate();
9863-
END REMOVE */
98649855
}
98659856

98669857
void Species::_ReadAncestralSequence(const char *p_file, Chromosome &p_chromosome)

0 commit comments

Comments
 (0)