From 49cdb4c7c8b3615f78bb39d972c523cd163231cd Mon Sep 17 00:00:00 2001 From: Jerome Kelleher Date: Thu, 1 Aug 2024 10:39:29 +0100 Subject: [PATCH] Fixup basic working C code lineage->tail --- lib/msprime.c | 53 ++++++++++++++++++++++++++------------- lib/tests/test_ancestry.c | 2 +- 2 files changed, 36 insertions(+), 19 deletions(-) diff --git a/lib/msprime.c b/lib/msprime.c index 085fe821e..229ada286 100644 --- a/lib/msprime.c +++ b/lib/msprime.c @@ -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; @@ -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; @@ -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); } @@ -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"); @@ -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; @@ -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; @@ -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 */ @@ -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( @@ -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) { @@ -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; @@ -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; @@ -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) { @@ -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; @@ -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); @@ -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; @@ -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; diff --git a/lib/tests/test_ancestry.c b/lib/tests/test_ancestry.c index 11833c1a7..82de438a7 100644 --- a/lib/tests/test_ancestry.c +++ b/lib/tests/test_ancestry.c @@ -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;