Skip to content

Commit

Permalink
Fixup basic working C code lineage->tail
Browse files Browse the repository at this point in the history
  • Loading branch information
@jeromekelleher
jeromekelleher committed Aug 1, 2024
1 parent c6b20da commit 49cdb4c
Show file tree
Hide file tree
Showing 2 changed files with 36 additions and 19 deletions.
53 changes: 35 additions & 18 deletions lib/msprime.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -892,8 +892,8 @@ lineage_reset_segments(lineage_t *self)
}

static lineage_t *MSP_WARN_UNUSED
msp_alloc_lineage(
msp_t *self, segment_t *head, segment_t *tail, population_id_t population, label_id_t label)
msp_alloc_lineage(msp_t *self, segment_t *head, segment_t *tail,
population_id_t population, label_id_t label)
{
lineage_t *lin = NULL;

Expand Down Expand Up @@ -1392,7 +1392,7 @@ msp_insert_hull(msp_t *self, hull_t *hull)
} else {
c = avl_search_closest(hulls_right, &query, &query_node);
/* query < node->item ==> c = -1 */
num_ending_before_left = (uint64_t) avl_index(query_node) + (uint64_t) (c != -1);
num_ending_before_left = (uint64_t) avl_index(query_node) + (uint64_t)(c != -1);
}
/* set number of pairs coalescing with hull */
count = num_starting_before_left - num_ending_before_left;
Expand Down Expand Up @@ -1895,7 +1895,7 @@ msp_verify_non_empty_populations(msp_t *self)

for (avl_node = self->non_empty_populations.head; avl_node != NULL;
avl_node = avl_node->next) {
j = (tsk_id_t) (intptr_t) avl_node->item;
j = (tsk_id_t)(intptr_t) avl_node->item;
tsk_bug_assert(msp_get_num_population_ancestors(self, j) > 0);
}

Expand Down Expand Up @@ -2280,7 +2280,7 @@ msp_print_state(msp_t *self, FILE *out)
}
fprintf(out, "non_empty_populations = [");
for (a = self->non_empty_populations.head; a != NULL; a = a->next) {
j = (uint32_t) (intptr_t) a->item;
j = (uint32_t)(intptr_t) a->item;
fprintf(out, "%d,", j);
}
fprintf(out, "]\n");
Expand Down Expand Up @@ -3347,7 +3347,8 @@ msp_recombination_event(msp_t *self, label_id_t label, segment_t **lhs, segment_
}
tsk_bug_assert(alpha->left < alpha->right);
msp_set_segment_mass(self, alpha);
right_lineage = msp_alloc_lineage(self, alpha, NULL, left_lineage->population, label);
right_lineage
= msp_alloc_lineage(self, alpha, NULL, left_lineage->population, label);
if (right_lineage == NULL) {
ret = MSP_ERR_NO_MEMORY;
goto out;
Expand Down Expand Up @@ -3424,7 +3425,7 @@ msp_gene_conversion_event(msp_t *self, label_id_t label)
{
int ret = 0;
segment_t *x, *y, *alpha, *head, *tail, *z, *new_individual_head;
lineage_t *new_lineage;
lineage_t *lineage, *new_lineage;
double left_breakpoint, right_breakpoint, tl;
bool insert_alpha;
hull_t *hull = NULL;
Expand All @@ -3441,7 +3442,8 @@ msp_gene_conversion_event(msp_t *self, label_id_t label)
goto out;
}

population = y->lineage->population;
lineage = y->lineage;
population = lineage->population;
x = y->prev;

/* generate tract length */
Expand Down Expand Up @@ -3611,12 +3613,17 @@ msp_gene_conversion_event(msp_t *self, label_id_t label)
new_individual_head = head;
}
if (new_individual_head != NULL) {
new_lineage = msp_alloc_lineage(self, new_individual_head, NULL, population, label);
new_lineage
= msp_alloc_lineage(self, new_individual_head, NULL, population, label);
if (new_lineage == NULL) {
ret = MSP_ERR_NO_MEMORY;
goto out;
}
lineage_reset_segments(new_lineage);
ret = msp_insert_individual(self, new_lineage);
if (ret != 0) {
goto out;
}
if (self->model.type == MSP_MODEL_SMC_K) {
tsk_bug_assert(tract_hull_left < tract_hull_right);
tract_hull_right = GSL_MIN(
Expand All @@ -3635,6 +3642,9 @@ msp_gene_conversion_event(msp_t *self, label_id_t label)
} else {
self->num_noneffective_gc_events++;
}

lineage_reset_segments(lineage);

if (self->additional_nodes & MSP_NODE_IS_GC_EVENT) {
ret = msp_store_arg_gene_conversion(self, tail, alpha, head);
if (ret != 0) {
Expand Down Expand Up @@ -4346,7 +4356,8 @@ msp_insert_root_segments(msp_t *self, const segment_t *head, segment_t **new_hea
}
copy->prev = prev;
if (prev == NULL) {
lineage = msp_alloc_lineage(self, copy, NULL, node_population[head->value], label);
lineage = msp_alloc_lineage(
self, copy, NULL, node_population[head->value], label);
if (lineage == NULL) {
ret = MSP_ERR_NO_MEMORY;
goto out;
Expand Down Expand Up @@ -5181,7 +5192,7 @@ msp_gene_conversion_left_event(msp_t *self, label_id_t label)
double h = gsl_rng_uniform(self->rng) * gc_left_total;
double tl, bp, lhs_old_right, lhs_new_right;
population_id_t population;
lineage_t *lineage;
lineage_t *lineage, *new_lineage;
segment_t *y, *x, *alpha;
hull_t *rhs_hull;
hull_t *lhs_hull = NULL;
Expand Down Expand Up @@ -5269,15 +5280,21 @@ msp_gene_conversion_left_event(msp_t *self, label_id_t label)
}
lhs_new_right = y->right;

lineage = msp_alloc_lineage(self, alpha, NULL, population, label);
if (lineage == NULL) {
new_lineage = msp_alloc_lineage(self, alpha, NULL, population, label);
if (new_lineage == NULL) {
ret = MSP_ERR_NO_MEMORY;
goto out;
}
msp_set_segment_mass(self, alpha);
tsk_bug_assert(alpha->prev == NULL);

ret = msp_insert_individual(self, lineage);
lineage_reset_segments(new_lineage);
lineage_reset_segments(lineage);

ret = msp_insert_individual(self, new_lineage);
if (ret != 0) {
goto out;
}
if (self->additional_nodes & MSP_NODE_IS_GC_EVENT) {
ret = msp_store_arg_gene_conversion(self, NULL, y, alpha);
if (ret != 0) {
Expand Down Expand Up @@ -5376,7 +5393,7 @@ msp_run_coalescent(msp_t *self, double max_time, unsigned long max_events)
ca_pop_id = 0;
for (avl_node = self->non_empty_populations.head; avl_node != NULL;
avl_node = avl_node->next) {
pop_id = (tsk_id_t) (intptr_t) avl_node->item;
pop_id = (tsk_id_t)(intptr_t) avl_node->item;
t_temp = self->get_common_ancestor_waiting_time(self, pop_id, label);
if (t_temp < ca_t_wait) {
ca_t_wait = t_temp;
Expand All @@ -5390,7 +5407,7 @@ msp_run_coalescent(msp_t *self, double max_time, unsigned long max_events)
mig_dest_pop = 0;
for (avl_node = self->non_empty_populations.head; avl_node != NULL;
avl_node = avl_node->next) {
pop_id_j = (tsk_id_t) (intptr_t) avl_node->item;
pop_id_j = (tsk_id_t)(intptr_t) avl_node->item;
pop = &self->populations[pop_id_j];
n = avl_count(&pop->ancestors[label]);
tsk_bug_assert(n > 0);
Expand Down Expand Up @@ -7505,7 +7522,7 @@ msp_instantaneous_bottleneck(msp_t *self, demographic_event_t *event)
for (u = 0; u < (tsk_id_t) n; u++) {
lineages[u] = u;
}
for (u = 0; u < (tsk_id_t) (2 * n); u++) {
for (u = 0; u < (tsk_id_t)(2 * n); u++) {
pi[u] = TSK_NULL;
}
j = 0;
Expand Down Expand Up @@ -8380,7 +8397,7 @@ genic_selection_generate_trajectory(sweep_t *self, msp_t *simulator,
alpha = 2 * pop_size * trajectory.s;
x = 1.0
- genic_selection_stochastic_forwards(
trajectory.dt, 1.0 - x, alpha, gsl_rng_uniform(rng));
trajectory.dt, 1.0 - x, alpha, gsl_rng_uniform(rng));
/* need our recored traj to stay in bounds */
t += trajectory.dt;
sim_time += trajectory.dt * pop_size * simulator->ploidy;
Expand Down
2 changes: 1 addition & 1 deletion lib/tests/test_ancestry.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1297,7 +1297,7 @@ test_mixed_hudson_smc(void)
/* Run for 1 event each for the models, interleaving */
j = 0;
while ((ret = msp_run(&msp, DBL_MAX, 1)) == MSP_EXIT_MAX_EVENTS) {
printf("X\n");
/* msp_print_state(&msp, stdout); */
msp_verify(&msp, 0);
CU_ASSERT_FALSE(msp_is_completed(&msp));
model = msp_get_model(&msp)->type;
Expand Down

0 comments on commit 49cdb4c

Please sign in to comment.