Skip to content

put binary derived state into the mutation table metadata #665

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
bhaller wants to merge 1 commit into
base: multitrait
Choose a base branch
from
+127 −103
Open

put binary derived state into the mutation table metadata #665

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 3 additions & 0 deletions VERSIONS
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -208,6 +208,9 @@ multitrait branch:
internally, split baselineOffset into a "from substitutions" component and a "from user" component, and expose those components in .trees metadata, to fix #661
the top-level JSON schema now has two keys, baselineOffsetFromUser and baselineOffsetFromSubstitutions, that provide these two components
SLiM writes both of them, for use on the Python side; when reading a .trees, SLiM uses baselineOffsetFromUser and ignores baselineOffsetFromSubstitutions (calculating that itself)
shift to putting binary derived state info in the mutation table metadata, in addition to the ASCII derived state info in its own column, when on disk (i.e., in .trees files); see #664
in-memory, the binary derived state info remains in the derived state column and the mutation table metadata is kept empty, as before
the schema for the mutation table metadata now encodes an array (of length determined by the metadata itself) of int64_t SLiM mutation IDs


version 5.2 (Eidos version 4.2):
Expand Down
34 changes: 33 additions & 1 deletion core/slim_globals.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2177,7 +2177,39 @@ const std::string gSLiM_tsk_metadata_binary_schema_FORMAT_SOURCE = R"V0G0N({

const std::string gSLiM_tsk_edge_metadata_schema_SOURCE = "";
const std::string gSLiM_tsk_site_metadata_schema_SOURCE = "";
const std::string gSLiM_tsk_mutation_metadata_schema_SOURCE = ""; // this is now managed by gSLiM_tsk_metadata_binary_schema in top-level metadata

// actual mutation metadata is now kept in the top-level metadata; see gSLiM_tsk_metadata_binary_schema
// the mutation metadata column now (on disk only) contains binary derived state info, in addition to the
// ASCII info kept in the derived state column; see DerivedStatesFromMetadata() and DerivedStatesToMetadata()
#pragma mark gSLiM_tsk_mutation_metadata_schema_SOURCE

const std::string gSLiM_tsk_mutation_metadata_schema_SOURCE = R"V0G0N({
"$schema": "http://json-schema.org/schema#",
"additionalProperties": false,
"codec": "struct",
"type": "object",
"description": "SLiM schema for representing binary derived state data in mutation metadata (the actual mutation metadata is stored in a table in top-level metadata now).",
"examples": [
{
"derived_states": [0, 1, 17]
}
],
"properties": {
"derived_states": {
"index": 1,
"type": "array",
"noLengthEncodingExhaustBuffer": true,
"description": "An array of SLiM mutation IDs (int64t), representing the (stacked) mutations contained by the derived state for the mutation. This information is also saved as a comma-separated ASCII string in the derived state column.",
"items": {
"binaryFormat": "q",
"type": "number"
}
}
},
"required": [
"derived_states"
]
})V0G0N";

// BCH 12/10/2024: Removed the type field, and changed the treatment of is_vacant. We have a tricky problem
// here, which is that is_vacant is now variable-length and there is no count. The number of byte (uint8_t)
Expand Down
189 changes: 89 additions & 100 deletions core/species.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8771,115 +8771,101 @@ void Species::CheckAutoSimplification(void)
}
}

void Species::DerivedStatesFromAscii(tsk_table_collection_t *p_tables)
void Species::DerivedStatesFromMetadata(tsk_table_collection_t *p_tables)
{
// This modifies p_tables in place, replacing the derived_state column of p_tables with a binary version.
tsk_mutation_table_t mutations_copy;
int ret = tsk_mutation_table_copy(&p_tables->mutations, &mutations_copy, 0);
if (ret < 0) handle_error("derived_from_ascii", ret);
// This is called when reading a .trees file. On disk, derived state information is kept in the mutation
// table's metadata column in binary, AND in the derived state column in ASCII (see #664 for discussion).
// Here we convert back to our in-memory format by swapping the derived state column's ASCII data with the
// metadata column's binary data, and then purging the ASCII data from the metadata column.

{
const char *derived_state = p_tables->mutations.derived_state;
tsk_size_t *derived_state_offset = p_tables->mutations.derived_state_offset;
std::vector<slim_mutationid_t> binary_derived_state;
std::vector<tsk_size_t> binary_derived_state_offset;
size_t derived_state_total_part_count = 0;

binary_derived_state_offset.emplace_back(0);

try {
for (size_t j = 0; j < p_tables->mutations.num_rows; j++)
{
std::string string_derived_state(derived_state + derived_state_offset[j], derived_state_offset[j+1] - derived_state_offset[j]);

if (string_derived_state.size() == 0)
{
// nothing to do for an empty derived state
}
else if (string_derived_state.find(',') == std::string::npos)
{
// a single mutation can be handled more efficiently, and this is the common case so it's worth optimizing
binary_derived_state.emplace_back((slim_mutationid_t)std::stoll(string_derived_state));
derived_state_total_part_count++;
}
else
{
// stacked mutations require that the derived state be separated to parse it
std::vector<std::string> derived_state_parts = Eidos_string_split(string_derived_state, ",");

for (std::string &derived_state_part : derived_state_parts)
binary_derived_state.emplace_back((slim_mutationid_t)std::stoll(derived_state_part));

derived_state_total_part_count += derived_state_parts.size();
}

binary_derived_state_offset.emplace_back((tsk_size_t)(derived_state_total_part_count * sizeof(slim_mutationid_t)));
}
} catch (...) {
EIDOS_TERMINATION << "ERROR (Species::DerivedStatesFromAscii): a mutation derived state was not convertible into an int64_t mutation id. The tree-sequence data may not be annotated for SLiM, or may be corrupted. If mutations were added in msprime, do you want to use the msprime.SLiMMutationModel?" << EidosTerminate();
}

if (binary_derived_state.size() == 0)
binary_derived_state.resize(1);

ret = tsk_mutation_table_set_columns(&p_tables->mutations,
mutations_copy.num_rows,
mutations_copy.site,
mutations_copy.node,
mutations_copy.parent,
mutations_copy.time,
(char *)binary_derived_state.data(),
binary_derived_state_offset.data(),
mutations_copy.metadata,
mutations_copy.metadata_offset);
if (ret < 0) handle_error("derived_from_ascii", ret);
}
// We want to do this efficiently, without making copies of buffers, etc., so we munge around in tskit's
// structs. Maybe there is a way to do this more cleanly and safely? :-O

assert(p_tables != nullptr);
assert(p_tables->mutations.derived_state != nullptr);
assert(p_tables->mutations.metadata != nullptr);
assert(p_tables->mutations.metadata_offset != nullptr);

// Swap the derived state column and the metadata column; after this, binary derived states will be in the
// derived state column and ASCII derived states will be in the metadata column.
std::swap(p_tables->mutations.derived_state_length, p_tables->mutations.metadata_length);
std::swap(p_tables->mutations.max_derived_state_length, p_tables->mutations.max_metadata_length);
std::swap(p_tables->mutations.max_derived_state_length_increment, p_tables->mutations.max_metadata_length_increment);
std::swap(p_tables->mutations.derived_state, p_tables->mutations.metadata);
std::swap(p_tables->mutations.derived_state_offset, p_tables->mutations.metadata_offset);

// Then empty out the metadata column, which we don't use while simulating; we want to free up the buffer
// for the metadata, rather than having the (perhaps large capacity) buffer hanging around forever. The
// only way to do that is by freeing the old buffer and mallocing a new minimal buffer.
free(p_tables->mutations.metadata);
p_tables->mutations.metadata = (char *)malloc(1); // avoids platform-dependencies on zero-length malloc
if (!p_tables->mutations.metadata)
EIDOS_TERMINATION << "ERROR (Species::DerivedStatesFromMetadata): allocation failed; you may need to raise the memory limit for SLiM." << EidosTerminate(nullptr);

tsk_mutation_table_free(&mutations_copy);
p_tables->mutations.metadata_length = 0;
p_tables->mutations.max_metadata_length = 1;

// Zero out all the metadata offsets; there is an entry for each row, plus one
EIDOS_BZERO(p_tables->mutations.metadata_offset, (p_tables->mutations.num_rows + 1) * sizeof(tsk_size_t));
}

void Species::DerivedStatesToAscii(tsk_table_collection_t *p_tables)
void Species::DerivedStatesToMetadata(tsk_table_collection_t *p_tables)
{
// This modifies p_tables in place, replacing the derived_state column of p_tables with an ASCII version.
tsk_mutation_table_t mutations_copy;
int ret = tsk_mutation_table_copy(&p_tables->mutations, &mutations_copy, 0);
if (ret < 0) handle_error("derived_to_ascii", ret);
// This is called when saving out a .trees file. On disk, derived state information is kept in the mutation
// table's metadata column in binary, AND in the derived state column in ASCII (see #664 for discussion).
// Here we convert from our in-memory format by swapping the derived state column's binary data with the
// empty metadata column, and then adding ASCII derived state information.

// We want to do this efficiently, without making copies of buffers, etc., so we munge around in tskit's
// structs. Maybe there is a way to do this more cleanly and safely? :-O

assert(p_tables != nullptr);
assert(p_tables->mutations.derived_state != nullptr);
assert(p_tables->mutations.metadata != nullptr);
assert(p_tables->mutations.metadata_offset != nullptr);

// Generate the ASCII derived state column, for use on disk, and put it into the metadata column.
const char *derived_state = p_tables->mutations.derived_state;
tsk_size_t *derived_state_offset = p_tables->mutations.derived_state_offset;
std::string text_derived_state;
tsk_size_t *text_derived_state_offset = p_tables->mutations.metadata_offset;

text_derived_state_offset[0] = 0;

for (size_t j = 0; j < p_tables->mutations.num_rows; j++)
{
const char *derived_state = p_tables->mutations.derived_state;
tsk_size_t *derived_state_offset = p_tables->mutations.derived_state_offset;
std::string text_derived_state;
std::vector<tsk_size_t> text_derived_state_offset;
slim_mutationid_t *int_derived_state = (slim_mutationid_t *)(derived_state + derived_state_offset[j]);
size_t cur_derived_state_length = (derived_state_offset[j+1] - derived_state_offset[j])/sizeof(slim_mutationid_t);

text_derived_state_offset.emplace_back(0);

for (size_t j = 0; j < p_tables->mutations.num_rows; j++)
for (size_t i = 0; i < cur_derived_state_length; i++)
{
slim_mutationid_t *int_derived_state = (slim_mutationid_t *)(derived_state + derived_state_offset[j]);
size_t cur_derived_state_length = (derived_state_offset[j+1] - derived_state_offset[j])/sizeof(slim_mutationid_t);

for (size_t i = 0; i < cur_derived_state_length; i++)
{
if (i != 0) text_derived_state.append(",");
text_derived_state.append(std::to_string(int_derived_state[i]));
}
text_derived_state_offset.emplace_back((tsk_size_t)text_derived_state.size());
if (i != 0) text_derived_state.append(",");
text_derived_state.append(std::to_string(int_derived_state[i]));
}

ret = tsk_mutation_table_set_columns(&p_tables->mutations,
mutations_copy.num_rows,
mutations_copy.site,
mutations_copy.node,
mutations_copy.parent,
mutations_copy.time,
text_derived_state.c_str(),
text_derived_state_offset.data(),
mutations_copy.metadata,
mutations_copy.metadata_offset);
if (ret < 0) handle_error("derived_to_ascii", ret);
text_derived_state_offset[j + 1] = (tsk_size_t)text_derived_state.size();
}

tsk_mutation_table_free(&mutations_copy);
// Copy the ASCII data into a new malloced block and replace any existing metadata with the ASCII.
// FIXME it'd be nice to do the work ourselves into a malloced buffer we own, to avoid the copy.
tsk_size_t metadata_size = text_derived_state.size() * sizeof(char);
char *new_metadata_buffer = (char *)malloc(metadata_size);
if (!new_metadata_buffer)
EIDOS_TERMINATION << "ERROR (Species::DerivedStatesToMetadata): allocation failed; you may need to raise the memory limit for SLiM." << EidosTerminate(nullptr);

memcpy(new_metadata_buffer, text_derived_state.c_str(), metadata_size);

free(p_tables->mutations.metadata);
p_tables->mutations.metadata = new_metadata_buffer;
p_tables->mutations.metadata_length = metadata_size;
p_tables->mutations.max_metadata_length = metadata_size;

// Swap the derived state column and the metadata column; after this, binary derived states will
// be in the metadata column, and ASCII derived states will be in the derived state column.
std::swap(p_tables->mutations.derived_state_length, p_tables->mutations.metadata_length);
std::swap(p_tables->mutations.max_derived_state_length, p_tables->mutations.max_metadata_length);
std::swap(p_tables->mutations.max_derived_state_length_increment, p_tables->mutations.max_metadata_length_increment);
std::swap(p_tables->mutations.derived_state, p_tables->mutations.metadata);
std::swap(p_tables->mutations.derived_state_offset, p_tables->mutations.metadata_offset);
}

void Species::AddIndividualsToTable(Individual * const *p_individual, size_t p_num_individuals, tsk_table_collection_t *p_tables, INDIVIDUALS_HASH *p_individuals_hash, tsk_flags_t p_flags)
Expand Down Expand Up @@ -10749,7 +10735,7 @@ void Species::WriteTreeSequence(std::string &p_recording_tree_path, bool p_simpl
// Write out the copied tables
{
// derived state data must be in ASCII (or unicode) on disk, according to tskit policy
DerivedStatesToAscii(&output_tables);
DerivedStatesToMetadata(&output_tables);

// In nucleotide-based models, put an ASCII representation of the reference sequence into the tables
if (nucleotide_based_)
Expand Down Expand Up @@ -13862,8 +13848,8 @@ slim_tick_t Species::_InitializePopulationFromTskitBinaryFile(const char *p_file

ReadTreeSequenceMetadata(treeSeqInfo, &metadata_tick, &metadata_cycle, &file_model_type, &file_version, mut_metadata_table);

// convert ASCII derived-state data, which is the required format on disk, back to our in-memory binary format
DerivedStatesFromAscii(&treeSeqInfo.tables_);
// shift derived state information from the mutation metadata column to the derived state column
DerivedStatesFromMetadata(&treeSeqInfo.tables_);

// in nucleotide-based models, read the ancestral sequence; we do this ourselves, directly from kastore, to avoid having
// tskit make a full ASCII copy of the reference sequences from kastore into tables_; see tsk_table_collection_load() above
Expand Down Expand Up @@ -14022,7 +14008,10 @@ slim_tick_t Species::_InitializePopulationFromTskitDirectory(std::string p_direc
EIDOS_TERMINATION << "ERROR (Species::_InitializePopulationFromTskitDirectory): the .trees files for chromosomes have different file versions (" << file_version << " versus " << this_file_version << "). This must be consistent across all files." << EidosTerminate();
}

DerivedStatesFromAscii(&treeSeqInfo.tables_);
// shift derived state information from the mutation metadata column to the derived state column
DerivedStatesFromMetadata(&treeSeqInfo.tables_);

// in nucleotide-based models, read the ancestral sequence
_ReadAncestralSequence(expected_path.c_str(), *chromosome);

// The first chromosome uses _InstantiateSLiMObjectsFromTables() and creates the subpopulations, etc.,
Expand Down
4 changes: 2 additions & 2 deletions core/species.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -728,8 +728,8 @@ class Species : public EidosDictionaryUnretained
void DisconnectCopiedSharedTables(tsk_table_collection_t &p_tables); // zeroes out the shared table copies in p_tables

static void handle_error(const std::string &msg, int error) __attribute__((__noreturn__)) __attribute__((cold)) __attribute__((analyzer_noreturn));
static void DerivedStatesFromAscii(tsk_table_collection_t *p_tables);
static void DerivedStatesToAscii(tsk_table_collection_t *p_tables);
static void DerivedStatesFromMetadata(tsk_table_collection_t *p_tables);
static void DerivedStatesToMetadata(tsk_table_collection_t *p_tables);

bool _SubpopulationIDInUse(slim_objectid_t p_subpop_id) const;
void RecordTablePosition(void);
Expand Down
Loading