@@ -118,30 +118,39 @@ tsk_haplotype_reset_bitset(const tsk_haplotype_t *self)
118118}
119119
120120static int
121- tsk_haplotype_build_child_index (tsk_haplotype_t * self )
121+ tsk_haplotype_build_parent_index (tsk_haplotype_t * self )
122122{
123123 int ret = 0 ;
124124 tsk_size_t j ;
125125 const tsk_table_collection_t * tables = self -> tree_sequence -> tables ;
126126 const tsk_edge_table_t * edges = & tables -> edges ;
127127 tsk_size_t num_edges = edges -> num_rows ;
128128 tsk_haplotype_edge_sort_t * sorted = NULL ;
129+ int32_t * offsets = NULL ;
129130
130131 if (num_edges == 0 ) {
131- self -> child_order = NULL ;
132- self -> child_offsets
133- = tsk_calloc (self -> num_nodes + 1 , sizeof (* self -> child_offsets ));
134- if (self -> child_offsets == NULL ) {
135- ret = tsk_trace_error (TSK_ERR_NO_MEMORY );
136- goto out ;
132+ self -> parent_edge_index = NULL ;
133+ if (self -> num_nodes > 0 ) {
134+ self -> parent_index_range
135+ = tsk_calloc (self -> num_nodes * 2 , sizeof (* self -> parent_index_range ));
136+ if (self -> parent_index_range == NULL ) {
137+ ret = tsk_trace_error (TSK_ERR_NO_MEMORY );
138+ goto out ;
139+ }
140+ } else {
141+ self -> parent_index_range = NULL ;
137142 }
138143 goto out ;
139144 }
140145
141146 sorted = tsk_malloc (num_edges * sizeof (* sorted ));
142- self -> child_order = tsk_malloc (num_edges * sizeof (* self -> child_order ));
143- self -> child_offsets = tsk_calloc (self -> num_nodes + 1 , sizeof (* self -> child_offsets ));
144- if (sorted == NULL || self -> child_order == NULL || self -> child_offsets == NULL ) {
147+ self -> parent_edge_index = tsk_malloc (num_edges * sizeof (* self -> parent_edge_index ));
148+ self -> parent_index_range
149+ = tsk_malloc (self -> num_nodes * 2 * sizeof (* self -> parent_index_range ));
150+ offsets = tsk_calloc (self -> num_nodes + 1 , sizeof (* offsets ));
151+ if (sorted == NULL || self -> parent_edge_index == NULL
152+ || (self -> num_nodes > 0 && self -> parent_index_range == NULL )
153+ || offsets == NULL ) {
145154 ret = tsk_trace_error (TSK_ERR_NO_MEMORY );
146155 goto out ;
147156 }
@@ -156,16 +165,26 @@ tsk_haplotype_build_child_index(tsk_haplotype_t *self)
156165 for (j = 0 ; j < num_edges ; j ++ ) {
157166 tsk_id_t child = sorted [j ].child ;
158167 if (child >= 0 && child < (tsk_id_t ) self -> num_nodes ) {
159- self -> child_offsets [child + 1 ]++ ;
168+ offsets [child + 1 ]++ ;
160169 }
161- self -> child_order [j ] = sorted [j ].edge_id ;
162170 }
163- for (j = 0 ; j < self -> num_nodes ; j ++ ) {
164- self -> child_offsets [j + 1 ] += self -> child_offsets [j ];
171+ for (j = 0 ; j < (tsk_size_t ) self -> num_nodes ; j ++ ) {
172+ offsets [j + 1 ] += offsets [j ];
173+ self -> parent_index_range [2 * j ] = offsets [j ];
174+ self -> parent_index_range [2 * j + 1 ] = offsets [j + 1 ];
175+ }
176+
177+ for (j = 0 ; j < num_edges ; j ++ ) {
178+ tsk_id_t child = sorted [j ].child ;
179+ if (child >= 0 && child < (tsk_id_t ) self -> num_nodes ) {
180+ int32_t pos = offsets [child ]++ ;
181+ self -> parent_edge_index [pos ] = sorted [j ].edge_id ;
182+ }
165183 }
166184
167185out :
168186 tsk_safe_free (sorted );
187+ tsk_safe_free (offsets );
169188 return ret ;
170189}
171190
@@ -427,7 +446,7 @@ tsk_haplotype_init(tsk_haplotype_t *self, const tsk_treeseq_t *tree_sequence,
427446 self -> num_edges = tables -> edges .num_rows ;
428447 self -> site_positions = sites -> position + site_start ;
429448
430- ret = tsk_haplotype_build_child_index (self );
449+ ret = tsk_haplotype_build_parent_index (self );
431450 if (ret != 0 ) {
432451 goto out ;
433452 }
@@ -523,10 +542,15 @@ tsk_haplotype_decode(tsk_haplotype_t *self, tsk_id_t node, int8_t *haplotype)
523542 }
524543 }
525544
526- int32_t child_start = self -> child_offsets [node ];
527- int32_t child_stop = self -> child_offsets [node + 1 ];
545+ int32_t child_start = 0 ;
546+ int32_t child_stop = 0 ;
547+ if (self -> parent_index_range != NULL ) {
548+ int32_t range_offset = node * 2 ;
549+ child_start = self -> parent_index_range [range_offset ];
550+ child_stop = self -> parent_index_range [range_offset + 1 ];
551+ }
528552 for (int32_t i = child_start ; i < child_stop ; i ++ ) {
529- tsk_id_t edge = self -> child_order [i ];
553+ tsk_id_t edge = self -> parent_edge_index [i ];
530554 int32_t start = self -> edge_start_index [edge ];
531555 int32_t end = self -> edge_end_index [edge ];
532556 if (start >= end ) {
@@ -565,11 +589,12 @@ tsk_haplotype_decode(tsk_haplotype_t *self, tsk_id_t node, int8_t *haplotype)
565589 }
566590
567591 parent_count = 0 ;
568- if (ancestor >= 0 ) {
569- child_start = self -> child_offsets [ancestor ];
570- child_stop = self -> child_offsets [ancestor + 1 ];
592+ if (ancestor >= 0 && self -> parent_index_range != NULL ) {
593+ int32_t range_offset = ancestor * 2 ;
594+ child_start = self -> parent_index_range [range_offset ];
595+ child_stop = self -> parent_index_range [range_offset + 1 ];
571596 for (int32_t i = child_start ; i < child_stop ; i ++ ) {
572- tsk_id_t parent_edge = self -> child_order [i ];
597+ tsk_id_t parent_edge = self -> parent_edge_index [i ];
573598 int32_t parent_start = self -> edge_start_index [parent_edge ];
574599 int32_t parent_end = self -> edge_end_index [parent_edge ];
575600 if (parent_start < interval_start ) {
@@ -593,6 +618,9 @@ tsk_haplotype_decode(tsk_haplotype_t *self, tsk_id_t node, int8_t *haplotype)
593618 parent_count ++ ;
594619 }
595620 }
621+ } else {
622+ child_start = 0 ;
623+ child_stop = 0 ;
596624 }
597625
598626 idx = tsk_haplotype_bitset_next (bits , self -> num_bit_words ,
@@ -638,8 +666,8 @@ tsk_haplotype_free(tsk_haplotype_t *self)
638666 tsk_safe_free (self -> node_mutation_offsets );
639667 tsk_safe_free (self -> node_mutation_sites );
640668 tsk_safe_free (self -> node_mutation_states );
641- tsk_safe_free (self -> child_order );
642- tsk_safe_free (self -> child_offsets );
669+ tsk_safe_free (self -> parent_edge_index );
670+ tsk_safe_free (self -> parent_index_range );
643671 tsk_safe_free (self -> edge_start_index );
644672 tsk_safe_free (self -> edge_end_index );
645673 tsk_safe_free (self -> edge_stack );
0 commit comments