diff --git a/Assembly.cpp b/Assembly.cpp index 0259cf3..407ee32 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -225,7 +225,6 @@ void* Build_hash_table(void* arg) int avalible_k = 0; - while (file_flag != 0) { @@ -802,13 +801,11 @@ overlap_region_alloc* overlap_list, int flag) } } - return available_overlaps; - } void get_new_candidates(long long readID, UC_Read* g_read, overlap_region_alloc* overlap_list, k_mer_pos_list_alloc* array_list, -HeapSq* heap, Candidates_list* l, double band_width_threshold, int keep_whole_chain) + Candidates_list* l, double band_width_threshold, int keep_whole_chain) { HPC_seq HPC_read; Hash_code k_code; @@ -819,7 +816,6 @@ HeapSq* heap, Candidates_list* l, double band_width_threshold, int keep_whole_ch uint64_t list_length; uint64_t sub_ID; - clear_Heap(heap); clear_Candidates_list(l); clear_k_mer_pos_list_alloc(array_list); @@ -827,18 +823,15 @@ HeapSq* heap, Candidates_list* l, double band_width_threshold, int keep_whole_ch recover_UC_Read(g_read, &R_INF, readID); - ///forward strand init_HPC_seq(&HPC_read, g_read->seq, g_read->length); init_Hash_code(&k_code); avalible_k = 0; - - while ((code = get_HPC_code(&HPC_read, &end_pos)) != 6) { if(code < 4) { - k_mer_append(&k_code,code, asm_opt.k_mer_length); + k_mer_append(&k_code, code, asm_opt.k_mer_length); avalible_k++; if (avalible_k >= asm_opt.k_mer_length) { @@ -857,38 +850,8 @@ HeapSq* heap, Candidates_list* l, double band_width_threshold, int keep_whole_ch } } + // BIG CHANGES WILL GO HERE!!! - - ///reverse complement strand - reverse_complement(g_read->seq, g_read->length); - init_HPC_seq(&HPC_read, g_read->seq, g_read->length); - init_Hash_code(&k_code); - avalible_k = 0; - - while ((code = get_HPC_code(&HPC_read, &end_pos)) != 6) - { - if(code < 4) - { - k_mer_append(&k_code,code, asm_opt.k_mer_length); - avalible_k++; - if (avalible_k >= asm_opt.k_mer_length) - { - list_length = locate_Total_Pos_Table(&PCB, &k_code, &list, asm_opt.k_mer_length, &sub_ID); - if (list_length != 0) - { - append_k_mer_pos_list_alloc(array_list, list, list_length, end_pos, 1); - } - } - } - else - { - avalible_k = 0; - init_Hash_code(&k_code); - } - } - - merge_k_mer_pos_list_alloc_heap_sort(array_list, l, heap); - calculate_overlap_region_by_chaining(l, overlap_list, readID, g_read->length, &R_INF, band_width_threshold, keep_whole_chain); } @@ -922,10 +885,6 @@ void* Overlap_calculate_heap_merge(void* arg) overlap_region_alloc overlap_list; init_overlap_region_alloc(&overlap_list); - HeapSq heap; - - Init_Heap(&heap); - Correct_dumy correct; init_Correct_dumy(&correct); @@ -946,7 +905,7 @@ void* Overlap_calculate_heap_merge(void* arg) for (i = thr_ID; i < (long long)R_INF.total_reads; i = i + asm_opt.thread_num) { ///get_new_candidates(i, &g_read, &overlap_list, &array_list, &heap, &l, THRESHOLD_RATE*1.5); - get_new_candidates(i, &g_read, &overlap_list, &array_list, &heap, &l, 0.02, 1); + get_new_candidates(i, &g_read, &overlap_list, &array_list, &l, 0.02, 1); clear_Cigar_record(¤t_cigar); clear_Round2_alignment(&second_round); @@ -982,7 +941,6 @@ void* Overlap_calculate_heap_merge(void* arg) destory_buffer_sub_block(¤t_sub_buffer); destory_Candidates_list(&l); destory_overlap_region_alloc(&overlap_list); - destory_Heap(&heap); destory_k_mer_pos_list_alloc(&array_list); destory_Graph(&POA_Graph); destory_Graph(&DAGCon); @@ -1040,9 +998,6 @@ void* Output_related_reads(void* arg) overlap_region_alloc overlap_list; init_overlap_region_alloc(&overlap_list); - HeapSq heap; - - Init_Heap(&heap); Correct_dumy correct; init_Correct_dumy(&correct); @@ -1072,7 +1027,7 @@ void* Output_related_reads(void* arg) memcmp(asm_opt.required_read_name, Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) { ////get_new_candidates(i, &g_read, &overlap_list, &array_list, &heap, &l, THRESHOLD_RATE*1.5); - get_new_candidates(i, &g_read, &overlap_list, &array_list, &heap, &l, 0.02, 1); + get_new_candidates(i, &g_read, &overlap_list, &array_list, &l, 0.02, 1); fprintf(stderr, ">%.*s\n", (int)Get_NAME_LENGTH((R_INF), i), Get_NAME((R_INF), i)); @@ -1097,7 +1052,6 @@ void* Output_related_reads(void* arg) destory_buffer_sub_block(¤t_sub_buffer); destory_Candidates_list(&l); destory_overlap_region_alloc(&overlap_list); - destory_Heap(&heap); destory_k_mer_pos_list_alloc(&array_list); destory_Graph(&POA_Graph); destory_Graph(&DAGCon); @@ -1626,9 +1580,6 @@ void* Final_overlap_calculate_heap_merge(void* arg) overlap_region_alloc overlap_list; init_overlap_region_alloc(&overlap_list); - HeapSq heap; - Init_Heap(&heap); - Cigar_record_alloc cigarline; init_Cigar_record_alloc(&cigarline); @@ -1636,15 +1587,11 @@ void* Final_overlap_calculate_heap_merge(void* arg) memset(c2n, 4, 256); c2n[(uint8_t)'A'] = c2n[(uint8_t)'a'] = 0; c2n[(uint8_t)'C'] = c2n[(uint8_t)'c'] = 1; c2n[(uint8_t)'G'] = c2n[(uint8_t)'g'] = 2; c2n[(uint8_t)'T'] = c2n[(uint8_t)'t'] = 3; // build the encoding table - - - - for (i = thr_ID; i < R_INF.total_reads; i = i + asm_opt.thread_num) { - get_new_candidates(i, &g_read, &overlap_list, &array_list, &heap, &l, 0.001, 0); + get_new_candidates(i, &g_read, &overlap_list, &array_list, &l, 0.001, 0); /** correct_overlap(&overlap_list, &R_INF, &g_read, &correct, &overlap_read, &POA_Graph, &DAGCon, &matched_overlap_0, &matched_overlap_1, &potiental_matched_overlap_0, &potiental_matched_overlap_1, @@ -1677,7 +1624,6 @@ void* Final_overlap_calculate_heap_merge(void* arg) destory_Candidates_list(&l); destory_overlap_region_alloc(&overlap_list); - destory_Heap(&heap); destory_k_mer_pos_list_alloc(&array_list); destory_UC_Read(&g_read); destory_UC_Read(&overlap_read); @@ -1946,15 +1892,3 @@ void Correct_Reads(int last_round) Correct_Reads(last_round - 1); } - - - - - - - - - - - - diff --git a/CommandLines.cpp b/CommandLines.cpp index 8734c47..ac424e6 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -63,7 +63,7 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->pat_index = NULL; asm_opt->mat_index = NULL; asm_opt->thread_num = 1; - asm_opt->k_mer_length = 40; + asm_opt->k_mer_length = 39; asm_opt->k_mer_min_freq = 3; asm_opt->k_mer_max_freq = 66; asm_opt->load_index_from_disk = 1; diff --git a/Hash_Table.cpp b/Hash_Table.cpp index 60d0bcf..9aa7f66 100644 --- a/Hash_Table.cpp +++ b/Hash_Table.cpp @@ -18,141 +18,6 @@ void overlap_region_sort_y_id(overlap_region *a, long long n) radix_sort_overlap_region_sort(a, a + n); } -void Init_Heap(HeapSq* HBT) -{ - HBT->MaxSize = 1000; - HBT->heap = (ElemType*)malloc(HBT->MaxSize*sizeof(ElemType)); - HBT->index_i = (uint64_t*)malloc(HBT->MaxSize*sizeof(uint64_t)); - HBT->len = 0; -} - -void destory_Heap(HeapSq* HBT) -{ - free(HBT->heap); - free(HBT->index_i); -} - -void clear_Heap(HeapSq* HBT) -{ - HBT->len = 0; -} - - - -inline int cmp_ElemType(ElemType* x, ElemType* y) -{ - - if (x->node.readID < y->node.readID) - { - return 1; - } - else if (x->node.readID > y->node.readID) - { - return 2; - } - else - { if (x->node.strand < y->node.strand) - { - return 1; - } - else if (x->node.strand > y->node.strand) - { - return 2; - } - else - { - - if(x->node.offset < y->node.offset) - { - return 1; - } - else if(x->node.offset > y->node.offset) - { - return 2; - } - else - { - if (x->node.self_offset < y->node.self_offset) - { - return 1; - } - else if (x->node.self_offset > y->node.self_offset) - { - return 2; - } - else ///if both r_pos_x and self_offset are equal, offset must be equal - { - return 0; - } - - } - } - } -} - - -inline void Insert_Heap(HeapSq* HBT, ElemType* x) -{ - long long i, j; - if (HBT->len == HBT->MaxSize) - { - HBT->MaxSize = 2*HBT->MaxSize; - HBT->heap = (ElemType*)realloc(HBT->heap, HBT->MaxSize*sizeof(ElemType)); - HBT->index_i = (uint64_t*)realloc(HBT->index_i,HBT->MaxSize*sizeof(uint64_t)); - } - HBT->heap[HBT->len] = *x; //add element to tail - HBT->len++; - i = HBT->len - 1; - while (i != 0) - { - j = (i - 1) / 2; - ///if (x >= HBT->heap[j]) - ///1: xheap[j])!=1) - break; - HBT->heap[i] = HBT->heap[j]; - i = j; - } - HBT->heap[i] = *x; -} - -inline int DeleteHeap(HeapSq* HBT, ElemType* get) -{ - ElemType temp, x; - int i, j; - if (HBT->len == 0) - { - return 0; - } - temp = HBT->heap[0]; - HBT->len--; - if (HBT->len == 0) - { - *get = temp; - return 2; - } - - x = HBT->heap[HBT->len]; - i = 0; - j = 2 * i + 1; - while (j <= HBT->len - 1) - { - ///if (j < HBT->len - 1 && HBT->heap[j] > HBT->heap[j+1]) - if (j < HBT->len - 1 && cmp_ElemType(&HBT->heap[j], &HBT->heap[j + 1]) == 2) - j++; - ///if (x <= HBT->heap[j]) - if (cmp_ElemType(&x, &HBT->heap[j])!=2) - break; - HBT->heap[i] = HBT->heap[j]; - i = j; - j = 2 * i + 1; - } - HBT->heap[i] = x; - - *get = temp; - return 1; -} - void init_overlap_region_alloc(overlap_region_alloc* list) { list->size = 1000; @@ -219,10 +84,8 @@ int get_fake_gap_shift(Fake_Cigar* x, int index) return result; } - - int append_inexact_overlap_region_alloc(overlap_region_alloc* list, overlap_region* tmp, -All_reads* R_INF, int add_beg_end) + All_reads* R_INF, int add_beg_end) { if (list->length + 1 > list->size) @@ -264,8 +127,6 @@ All_reads* R_INF, int add_beg_end) long long x_right_length = Get_READ_LENGTH((*R_INF), tmp->x_id) - tmp->x_pos_e - 1; long long y_right_length = Get_READ_LENGTH((*R_INF), tmp->y_id) - tmp->y_pos_e - 1; - - if(x_right_length <= y_right_length) { tmp->x_pos_e = Get_READ_LENGTH((*R_INF), tmp->x_id) - 1; @@ -276,7 +137,6 @@ All_reads* R_INF, int add_beg_end) tmp->x_pos_e = tmp->x_pos_e + y_right_length; tmp->y_pos_e = Get_READ_LENGTH((*R_INF), tmp->y_id) - 1; } - if (tmp->x_pos_strand == 1) { @@ -290,9 +150,6 @@ All_reads* R_INF, int add_beg_end) list->list[list->length].y_pos_s = Get_READ_LENGTH((*R_INF), tmp->y_id) - tmp->y_pos_e - 1; list->list[list->length].y_pos_strand = 1; - - - resize_fake_cigar(&(list->list[list->length].f_cigar), (tmp->f_cigar.length + 2)); if(add_beg_end == 1) { @@ -388,7 +245,6 @@ All_reads* R_INF, int add_beg_end) return 1; } - void append_overlap_region_alloc_debug(overlap_region_alloc* list, overlap_region* tmp) { @@ -414,7 +270,6 @@ void append_overlap_region_alloc_debug(overlap_region_alloc* list, overlap_regio list->length++; } - int cmp_by_x_pos_s(const void * a, const void * b) { if ((*(overlap_region*)a).x_pos_s > (*(overlap_region*)b).x_pos_s) @@ -474,7 +329,6 @@ int cmp_by_x_pos_e(const void * a, const void * b) } } - void debug_chain(k_mer_hit* a, long long a_n, Chain_Data* dp) { long long i, j, current_j; @@ -493,7 +347,7 @@ void debug_chain(k_mer_hit* a, long long a_n, Chain_Data* dp) if(j != -1) { distance_self_pos = a[current_j].self_offset - a[j].self_offset; - distance_pos = a[current_j].offset - a[j].offset; + distance_pos = ha_hit_get_offset(&a[current_j]) - ha_hit_get_offset(&a[j]); distance_gap = distance_pos > distance_self_pos? distance_pos - distance_self_pos : distance_self_pos - distance_pos; indels += distance_gap; @@ -517,7 +371,7 @@ void debug_chain(k_mer_hit* a, long long a_n, Chain_Data* dp) } long long get_chainLen(long long x_beg, long long x_end, long long xLen, -long long y_beg, long long y_end, long long yLen) + long long y_beg, long long y_end, long long yLen) { if(x_beg <= y_beg) { @@ -551,7 +405,7 @@ long long y_beg, long long y_end, long long yLen) ///double band_width_threshold = 0.05; void chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* result, -double band_width_threshold, int max_skip, int x_readLen, int y_readLen) + double band_width_threshold, int max_skip, int x_readLen, int y_readLen) { long long i, j; long long self_pos, pos, max_j, max_i, max_score, score, n_skip; @@ -567,7 +421,7 @@ double band_width_threshold, int max_skip, int x_readLen, int y_readLen) // fill the score and backtrack arrays for (i = 0; i < a_n; ++i) { - pos = a[i].offset; + pos = ha_hit_get_offset(&a[i]); self_pos = a[i].self_offset; max_j = -1; max_score = min_score; @@ -579,7 +433,7 @@ double band_width_threshold, int max_skip, int x_readLen, int y_readLen) ///may have a pre-cut condition for j for (j = i - 1; j >= 0; --j) { - distance_pos = pos - a[j].offset; + distance_pos = pos - ha_hit_get_offset(&a[j]); distance_self_pos = self_pos - a[j].self_offset; ///a has been sorted by a[].offset ///note for a, we do not have any two elements that have both equal offsets and self_offsets @@ -636,12 +490,8 @@ double band_width_threshold, int max_skip, int x_readLen, int y_readLen) dp->self_length[i] = max_self_length; } - ///debug_chain(a, a_n, dp); - - - max_score = -1; max_i = -1; long long mini_xLen = x_readLen * 2 + 2, tmp_xLen; @@ -652,12 +502,12 @@ double band_width_threshold, int max_skip, int x_readLen, int y_readLen) max_score = dp->score[i]; max_i = i; mini_xLen = get_chainLen(a[i].self_offset, a[i].self_offset, x_readLen, - a[i].offset, a[i].offset, y_readLen); + ha_hit_get_offset(&a[i]), ha_hit_get_offset(&a[i]), y_readLen); } else if(dp->score[i] == max_score) { tmp_xLen = get_chainLen(a[i].self_offset, a[i].self_offset, x_readLen, - a[i].offset, a[i].offset, y_readLen); + ha_hit_get_offset(&a[i]), ha_hit_get_offset(&a[i]), y_readLen); if(tmp_xLen < mini_xLen) { @@ -670,17 +520,16 @@ double band_width_threshold, int max_skip, int x_readLen, int y_readLen) } - clear_fake_cigar(&(result->f_cigar)); ///not a has been sorted by offset, that means has been sorted by query offset i = max_i; result->x_pos_e = a[i].self_offset; - result->y_pos_e = a[i].offset; + result->y_pos_e = ha_hit_get_offset(&a[i]); result->shared_seed = max_score; result->overlapLen = mini_xLen; distance_self_pos = result->x_pos_e - a[i].self_offset; - distance_pos = result->y_pos_e - a[i].offset; + distance_pos = result->y_pos_e - ha_hit_get_offset(&a[i]); long long pre_distance_gap = distance_pos - distance_self_pos; ///record first site ///the length of f_cigar should be at least 1 @@ -692,7 +541,7 @@ double band_width_threshold, int max_skip, int x_readLen, int y_readLen) while (i >= 0) { distance_self_pos = result->x_pos_e - a[i].self_offset; - distance_pos = result->y_pos_e - a[i].offset; + distance_pos = result->y_pos_e - ha_hit_get_offset(&a[i]); distance_gap = distance_pos - distance_self_pos; if(distance_gap != pre_distance_gap) { @@ -703,7 +552,7 @@ double band_width_threshold, int max_skip, int x_readLen, int y_readLen) chainLen++; result->x_pos_s = a[i].self_offset; - result->y_pos_s = a[i].offset; + result->y_pos_s = ha_hit_get_offset(&a[i]); i = dp->pre[i]; } } @@ -713,7 +562,7 @@ double band_width_threshold, int max_skip, int x_readLen, int y_readLen) while (i >= 0) { distance_self_pos = result->x_pos_e - a[i].self_offset; - distance_pos = result->y_pos_e - a[i].offset; + distance_pos = result->y_pos_e - ha_hit_get_offset(&a[i]); distance_gap = distance_pos - distance_self_pos; if(distance_gap == pre_distance_gap) { @@ -728,16 +577,14 @@ double band_width_threshold, int max_skip, int x_readLen, int y_readLen) chainLen++; result->x_pos_s = a[i].self_offset; - result->y_pos_s = a[i].offset; + result->y_pos_s = ha_hit_get_offset(&a[i]); i = dp->pre[i]; } } } - - void calculate_overlap_region_by_chaining(Candidates_list* candidates, overlap_region_alloc* overlap_list, -uint64_t readID, uint64_t readLength, All_reads* R_INF, double band_width_threshold, int add_beg_end) + uint64_t readID, uint64_t readLength, All_reads* R_INF, double band_width_threshold, int add_beg_end) { overlap_region tmp_region; long long i = 0; @@ -757,8 +604,8 @@ uint64_t readID, uint64_t readLength, All_reads* R_INF, double band_width_thresh i = 0; while (i < candidates->length) { - current_ID = candidates->list[i].readID; - current_stand = candidates->list[i].strand; + current_ID = ha_hit_get_readID(&candidates->list[i]); + current_stand = ha_hit_get_rev(&candidates->list[i]); ///reference read tmp_region.x_id = readID; @@ -776,9 +623,9 @@ uint64_t readID, uint64_t readLength, All_reads* R_INF, double band_width_thresh while (i < candidates->length && - current_ID == candidates->list[i].readID + current_ID == ha_hit_get_readID(&candidates->list[i]) && - current_stand == candidates->list[i].strand) + current_stand == ha_hit_get_rev(&candidates->list[i])) { sub_region_end = i; i++; @@ -811,10 +658,8 @@ uint64_t readID, uint64_t readLength, All_reads* R_INF, double band_width_thresh qsort(overlap_list->list, overlap_list->length, sizeof(overlap_region), cmp_by_x_pos_s); } - - void append_window_list(overlap_region* region, uint64_t x_start, uint64_t x_end, int y_start, int y_end, int error, -int extra_begin, int extra_end, int error_threshold) + int extra_begin, int extra_end, int error_threshold) { long long length = region->x_pos_e - region->x_pos_s + 1; @@ -843,8 +688,6 @@ int extra_begin, int extra_end, int error_threshold) region->w_list_length++; } - - void init_k_mer_pos_list_alloc(k_mer_pos_list_alloc* list) { list->size = 1000; @@ -864,9 +707,8 @@ void destory_k_mer_pos_list_alloc(k_mer_pos_list_alloc* list) } void append_k_mer_pos_list_alloc(k_mer_pos_list_alloc* list, k_mer_pos* n_list, uint64_t n_length, -uint64_t n_end_pos, uint8_t n_direction) + uint64_t n_end_pos, uint8_t n_direction) { - if (list->length + 1 > list->size) { list->size = list->size * 2; @@ -881,9 +723,6 @@ uint64_t n_end_pos, uint8_t n_direction) list->length++; } - - - int cmp_k_mer_pos_list(const void * a, const void * b) { if ((*(k_mer_pos_list*)a).length > (*(k_mer_pos_list*)b).length) @@ -898,18 +737,8 @@ int cmp_k_mer_pos_list(const void * a, const void * b) { return 0; } - } -inline void append_pos_to_Candidates_list(Candidates_list* candidates, ElemType* x) -{ - candidates->list[candidates->length] = x->node; - candidates->length++; -} - - - - void test_single_list(Candidates_list* candidates, k_mer_pos* n_list, uint64_t n_lengh, uint64_t end_pos, uint64_t strand) { uint64_t i; @@ -920,13 +749,13 @@ void test_single_list(Candidates_list* candidates, k_mer_pos* n_list, uint64_t n for (; j < candidates->length; j++) { if ( - n_list[i].offset == (uint64_t)candidates->list[j].offset + n_list[i].offset == (uint64_t)ha_hit_get_offset(&candidates->list[j]) && - n_list[i].readID == candidates->list[j].readID + n_list[i].readID == ha_hit_get_readID(&candidates->list[j]) && end_pos == (uint64_t)candidates->list[j].self_offset && - strand == candidates->list[j].strand + strand == ha_hit_get_rev(&candidates->list[j]) ) { break; @@ -938,82 +767,8 @@ void test_single_list(Candidates_list* candidates, k_mer_pos* n_list, uint64_t n fprintf(stderr, "ERROR 4\n"); } } - } - -void merge_k_mer_pos_list_alloc_heap_sort(k_mer_pos_list_alloc* list, Candidates_list* candidates, HeapSq* HBT) -{ - clear_Heap(HBT); - - uint64_t total_length = 0; - uint64_t i; - ElemType x, y; - //add the first element of each list to stack - for (i = 0; i < list->length; i++) - { - x.ID = i; - x.node.offset = list->list[i].list[0].offset; - x.node.readID = list->list[i].list[0].readID; - x.node.self_offset = list->list[i].end_pos; - x.node.strand = list->list[i].direction; - - Insert_Heap(HBT, &x); - - HBT->index_i[i] = 1; - - total_length = total_length + list->list[i].length; - } - - - candidates->length = 0; - if(total_length > (uint64_t)candidates->size) - { - candidates->size = total_length; - candidates->list = (k_mer_hit*)realloc(candidates->list, sizeof(k_mer_hit)*candidates->size); - candidates->tmp = (k_mer_hit*)realloc(candidates->tmp, sizeof(k_mer_hit)*candidates->size); - } - - uint64_t ID; - int flag; - ///while (flag = DeleteHeap(HBT, &x)) - while ((flag = DeleteHeap(HBT, &x))) - { - append_pos_to_Candidates_list(candidates, &x); - - i = HBT->index_i[x.ID]; - ID = x.ID; - - if (flag == 2) - { - for (; i < list->list[x.ID].length; i++) - { - y.ID = ID; - y.node.offset = list->list[x.ID].list[i].offset; - y.node.readID = list->list[x.ID].list[i].readID; - y.node.self_offset = list->list[x.ID].end_pos; - y.node.strand = list->list[x.ID].direction; - append_pos_to_Candidates_list(candidates, &y); - } - - break; - } - - if (i < list->list[x.ID].length) - { - y.ID = ID; - y.node.offset = list->list[x.ID].list[i].offset; - y.node.readID = list->list[x.ID].list[i].readID; - y.node.self_offset = list->list[x.ID].end_pos; - y.node.strand = list->list[x.ID].direction; - Insert_Heap(HBT, &y); - HBT->index_i[ID]++; - } - } -} - - - void init_Count_Table(Count_Table** table) { *table = ha_ct_init(); @@ -1233,127 +988,6 @@ void insert_H_peaks(H_peaks* h, long long index, long long value) h->list[index] += value; } -inline void RC_Hash_code(Hash_code* code, Hash_code* rc_code, int k) -{ - rc_code->x[0] = 0; - rc_code->x[1] = 0; - int i; - for (i = 0; i < k; i++) - { - rc_code->x[0] = rc_code->x[0] << 1; - rc_code->x[1] = rc_code->x[1] << 1; - rc_code->x[0] |= (((uint64_t)((code->x[0] >> i) & 1))^((uint64_t)1)); - rc_code->x[1] |= (((uint64_t)((code->x[1] >> i) & 1))^((uint64_t)1)); - } -} - - - -void get_peak_debug(Total_Count_Table* TCB, long long* min, long long* max) -{ - int i; - Count_Table* h; - khint_t k; - long long c_count; - H_peaks LH; - LH.list = NULL; - LH.length = 0; - uint64_t sub_ID; - uint64_t sub_key; - Hash_code code, rc_code, debug_code; - char str[100]; - char rc_str[100]; - - - - for (i = 0; i < TCB->size; i++) - { - h = TCB->sub_h[i]; - for (k = 0; k != kh_end(h); ++k) - { - if (kh_exist(h, k)) // test if a bucket contains data - { - sub_ID = i; - sub_key = kh_key(h, k); - - recover_hash_code(sub_ID, sub_key, &code, TCB->suffix_mode, - TCB->suffix_bits, asm_opt.k_mer_length); - RC_Hash_code(&code, &rc_code, asm_opt.k_mer_length); - RC_Hash_code(&rc_code, &debug_code, asm_opt.k_mer_length); - if(code.x[0] != debug_code.x[0] || code.x[1] != debug_code.x[1]) - { - fprintf(stderr, "error\n"); - } - - Hashcode_to_string(&code, str, asm_opt.k_mer_length); - Hashcode_to_string(&rc_code, rc_str, asm_opt.k_mer_length); - reverse_complement(str, asm_opt.k_mer_length); - if(memcmp(str, rc_str, asm_opt.k_mer_length) != 0) - { - fprintf(stderr, "error\n"); - int j; - for (j = 0; j < asm_opt.k_mer_length; j++) - { - fprintf(stderr, "%c",str[j]); - } - fprintf(stderr, "\n"); - - for (j = 0; j < asm_opt.k_mer_length; j++) - { - fprintf(stderr, "%c",rc_str[j]); - } - fprintf(stderr, "\n"); - - } - - - ///get_Total_Count_Table(&TCB, &k_code, k_mer_length); - - c_count = kh_val(h, k); - - if(get_Total_Count_Table(TCB, &code, asm_opt.k_mer_length) != c_count) - { - fprintf(stderr, "error\n"); - } - - insert_H_peaks(&LH, c_count, c_count); - } - } - } - - (*max) = -1; - (*min) = -1; - long long max_value = -1; - for (i = 0; i < (long long)LH.length; i++) - { - if(LH.list[i] >= max_value) - { - max_value = LH.list[i]; - (*max) = i; - } - } - - long long min_value = max_value; - for (i = 0; i < (long long)LH.length; i++) - { - if(LH.list[i] < min_value && LH.list[i] != 0) - { - min_value = LH.list[i]; - (*min) = i; - } - } - - for (i = 0; i < (long long)LH.length; i++) - { - ///fprintf(stderr, "%d, %d\n", i, LH.list[i]); - fprintf(stderr, "%lld\n", LH.list[i]); - } - - - - free(LH.list); -} - ///1: a > b; -1: a < b; 0: a=b int cmp_Hash_code(Hash_code* a, Hash_code* b) { @@ -1378,46 +1012,10 @@ int cmp_Hash_code(Hash_code* a, Hash_code* b) return 0; } -int get_total_freq(Total_Count_Table* TCB, uint64_t sub_ID, uint64_t sub_key, long long* T_count) +void get_total_freq(Total_Count_Table* TCB, uint64_t sub_ID, uint64_t sub_key, long long* T_count) { - Hash_code code, rc_code; - long long count, rc_count; - - recover_hash_code(sub_ID, sub_key, &code, TCB->suffix_mode, - TCB->suffix_bits, asm_opt.k_mer_length); - RC_Hash_code(&code, &rc_code, asm_opt.k_mer_length); - - count = get_Total_Count_Table(TCB, &code, asm_opt.k_mer_length); - rc_count = get_Total_Count_Table(TCB, &rc_code, asm_opt.k_mer_length); - (*T_count) = count + rc_count; - - if(count == 0) - { - return 0; - }///count > 0 && rc_count == 0 - else if(rc_count == 0) - { - return 1; - }///count > 0 && rc_count > 0 - else - { - int flag = cmp_Hash_code(&code, &rc_code); - - ///code > rc_code - if(flag > 0) - { - return 1; - }///code < rc_code - else if(flag < 0) - { - return 0; - } - else - { - (*T_count) = (*T_count)/2; - return 1; - } - } + khint_t t = ha_ct_get(TCB->sub_h[sub_ID], sub_key); + *T_count = kh_val(TCB->sub_h[sub_ID], t); } void get_peak(Total_Count_Table* TCB, long long* min, long long* max, long long* up_boundary) @@ -1432,8 +1030,6 @@ void get_peak(Total_Count_Table* TCB, long long* min, long long* max, long long* uint64_t sub_ID; uint64_t sub_key; - - for (i = 0; i < TCB->size; i++) { h = TCB->sub_h[i]; @@ -1444,10 +1040,8 @@ void get_peak(Total_Count_Table* TCB, long long* min, long long* max, long long* sub_ID = i; sub_key = kh_key(h, k); - if(get_total_freq(TCB, sub_ID, sub_key, &count)==1) - { - insert_H_peaks(&LH, count, count); - } + get_total_freq(TCB, sub_ID, sub_key, &count); + insert_H_peaks(&LH, count, count); } } } @@ -1492,9 +1086,6 @@ void get_peak(Total_Count_Table* TCB, long long* min, long long* max, long long* (*up_boundary) = (*max) * 10; } - - - long long min_value = max_value; //// seed with freq 1 is useless for (i = 2; i < (long long)LH.length && i < (*max); i++) @@ -1509,8 +1100,6 @@ void get_peak(Total_Count_Table* TCB, long long* min, long long* max, long long* free(LH.list); } - - void Traverse_Counting_Table(Total_Count_Table* TCB, Total_Pos_Table* PCB, int k_mer_min_freq, int k_mer_max_freq) { int i; @@ -1571,7 +1160,6 @@ void Traverse_Counting_Table(Total_Count_Table* TCB, Total_Pos_Table* PCB, int k ///kh_val(PCB->sub_h[sub_ID], t)++; fprintf(stderr, "ERROR\n"); } - PCB->useful_k_mer++; PCB->total_occ = PCB->total_occ + kh_val(h, k); @@ -1580,7 +1168,6 @@ void Traverse_Counting_Table(Total_Count_Table* TCB, Total_Pos_Table* PCB, int k } } - // fprintf(stdout, "useful_k_mer: %lld\n",PCB->useful_k_mer); // fprintf(stdout, "total_occ: %lld\n",PCB->total_occ); @@ -1615,13 +1202,8 @@ void Traverse_Counting_Table(Total_Count_Table* TCB, Total_Pos_Table* PCB, int k PCB->pos = (k_mer_pos*)malloc(sizeof(k_mer_pos)*PCB->total_occ); memset(PCB->pos, 0, sizeof(k_mer_pos)*PCB->total_occ); - - } - - - int cmp_k_mer_pos(const void * a, const void * b) { if ((*(k_mer_pos*)a).readID != (*(k_mer_pos*)b).readID) @@ -1657,7 +1239,6 @@ void clear_Chain_Data(Chain_Data* x) x->length = 0; } - void destory_Chain_Data(Chain_Data* x) { free(x->score); @@ -1666,7 +1247,6 @@ void destory_Chain_Data(Chain_Data* x) free(x->self_length); } - void resize_Chain_Data(Chain_Data* x, long long size) { if(size > x->size) @@ -1679,23 +1259,18 @@ void resize_Chain_Data(Chain_Data* x, long long size) } } - void init_Candidates_list(Candidates_list* l) { l->length = 0; l->size = 0; l->list = NULL; l->tmp = NULL; - l->foward_pos = 0; - l->rc_pos = 0; init_Chain_Data(&(l->chainDP)); } void clear_Candidates_list(Candidates_list* l) { l->length = 0; - l->foward_pos = 0; - l->rc_pos = 0; clear_Chain_Data(&(l->chainDP)); } @@ -1706,52 +1281,6 @@ void destory_Candidates_list(Candidates_list* l) destory_Chain_Data(&(l->chainDP)); } - - -int cmp_candidates_list(const void * a, const void * b) -{ - long long r_pos_a, r_pos_b; - - if((*((k_mer_hit*)a)).strand != (*((k_mer_hit*)b)).strand) - { - return (*((k_mer_hit*)a)).strand > (*((k_mer_hit*)b)).strand ? 1: -1; - } - else - { - if((*((k_mer_hit*)a)).readID != (*((k_mer_hit*)b)).readID) - { - return (*((k_mer_hit*)a)).readID > (*((k_mer_hit*)b)).readID ? 1: -1; - } - else - { - r_pos_a = (*((k_mer_hit*)a)).offset - (*((k_mer_hit*)a)).self_offset; - r_pos_b = (*((k_mer_hit*)b)).offset - (*((k_mer_hit*)b)).self_offset; - - if(r_pos_a != r_pos_b) - { - return r_pos_a > r_pos_b ? 1: -1; - } - else - { - if((*((k_mer_hit*)a)).self_offset != (*((k_mer_hit*)b)).self_offset) - { - return (*((k_mer_hit*)a)).self_offset > (*((k_mer_hit*)b)).self_offset ? 1: -1; - } - else - { - return 0; - } - - } - - } - - } -} - - - - void init_fake_cigar(Fake_Cigar* x) { x->buffer = NULL; @@ -1847,7 +1376,5 @@ void resize_window_list_alloc(window_list_alloc* x, long long size) { x->buffer[i].error = -1; } - - x->length = 0; } diff --git a/Hash_Table.h b/Hash_Table.h index a97ad81..458312d 100644 --- a/Hash_Table.h +++ b/Hash_Table.h @@ -35,8 +35,7 @@ KHASHL_MAP_INIT(static inline, Pos_Table, ha_pt, uint64_t, uint64_t, kh_hash_dum typedef struct { volatile int lock; - -}Hash_table_spin_lock; +} Hash_table_spin_lock; typedef struct { @@ -72,12 +71,11 @@ typedef struct uint64_t length; } k_mer_pos_list_alloc; - typedef struct { - int C_L[CIGAR_MAX_LENGTH]; - char C_C[CIGAR_MAX_LENGTH]; - int length; + int C_L[CIGAR_MAX_LENGTH]; + char C_C[CIGAR_MAX_LENGTH]; + int length; } CIGAR; typedef struct @@ -100,15 +98,14 @@ typedef struct window_list* buffer; long long length; long long size; -}window_list_alloc; - +} window_list_alloc; typedef struct { uint64_t* buffer; uint64_t length; uint64_t size; -}Fake_Cigar; +} Fake_Cigar; typedef struct { @@ -140,7 +137,6 @@ typedef struct window_list_alloc boundary_cigars; } overlap_region; - typedef struct { overlap_region* list; @@ -152,29 +148,24 @@ typedef struct typedef struct { - ///uint64_t offset; - long long offset; - ///uint64_t self_offset; - long long self_offset; - uint64_t readID; - uint8_t strand; + uint64_t opos; + uint32_t self_offset; // offset on the target read } k_mer_hit; - -typedef struct +static inline uint32_t ha_hit_get_readID(const k_mer_hit *h) { - k_mer_hit node; - uint64_t ID; -} ElemType; + return h->opos >> 33; +} - -typedef struct +static inline uint32_t ha_hit_get_rev(const k_mer_hit *h) { - ElemType* heap; - uint64_t* index_i; - int len; - int MaxSize; -} HeapSq; + return h->opos >> 32 & 1; +} + +static inline uint32_t ha_hit_get_offset(const k_mer_hit *h) +{ + return (uint32_t)h->opos; +} typedef struct { @@ -192,8 +183,6 @@ typedef struct k_mer_hit* tmp; long long length; long long size; - uint64_t foward_pos; - uint64_t rc_pos; Chain_Data chainDP; } Candidates_list; @@ -212,24 +201,9 @@ typedef struct uint64_t* k_mer_index; } Total_Pos_Table; - - - -inline uint64_t mod_d(uint64_t h_key, uint64_t low_key, uint64_t d) -{ - uint64_t result = (h_key >> 32) % d; - result = ((result << 32) + (h_key & (uint64_t)0xffffffff)) % d; - result = ((result << 32) + (low_key >> 32)) % d; - result = ((result << 32) + (low_key & (uint64_t)0xffffffff)) % d; - - return result; -} - - - ////suffix_bits = 64 in default inline int recover_hash_code(uint64_t sub_ID, uint64_t sub_key, Hash_code* code, -uint64_t suffix_mode, int suffix_bits, int k) + uint64_t suffix_mode, int suffix_bits, int k) // FIXME: not working right now { uint64_t h_key, low_key; h_key = low_key = 0; @@ -249,31 +223,39 @@ uint64_t suffix_mode, int suffix_bits, int k) return 1; } -///inline int get_sub_table(uint64_t* get_sub_ID, uint64_t* get_sub_key, Total_Count_Table* TCB, Hash_code* code, int k) -inline int get_sub_table(uint64_t* get_sub_ID, uint64_t* get_sub_key, uint64_t suffix_mode, int suffix_bits, -Hash_code* code, int k) +static inline int ha_code2rev(const Hash_code *code) { - uint64_t h_key, low_key; - ///k might be 64，so it is unsafe - ///low_key = code->x[0] | (code->x[1] << k); - low_key = code->x[0] | (code->x[1] << SAFE_SHIFT(k)); - //k cannot be 0, so this shift is safe - h_key = code->x[1] >> (64 - k); - - if(mod_d(h_key, low_key, MODE_VALUE) > 3) - { - return 0; - } - - uint64_t sub_ID = (low_key >> SAFE_SHIFT(suffix_bits)) | (h_key << (64 - suffix_bits)); - uint64_t sub_key = (low_key & suffix_mode); - - *get_sub_ID = sub_ID; - *get_sub_key = sub_key; - - return 1; + return code->x[1] < code->x[3]? 0 : 1; } +static inline int ha_get_sub_table_short(uint64_t* get_sub_ID, uint64_t* get_sub_key, uint64_t suffix_mode, int suffix_bits, Hash_code* code, int k) +{ // for k < 32 + int j = ha_code2rev(code); + uint64_t y = code->x[j<<1|1] << k | code->x[j<<1|0]; + y = yak_hash64(y, (1ULL<<(k+k)) - 1); + if (y % MODE_VALUE > 3) return 0; + *get_sub_ID = y >> suffix_bits; + *get_sub_key = y & suffix_mode; + return 1; +} + +static inline int ha_get_sub_table_long(uint64_t* get_sub_ID, uint64_t* get_sub_key, uint64_t suffix_mode, int suffix_bits, Hash_code* code, int k) +{ // for k > 32 + int j = ha_code2rev(code); + uint64_t y = code->x[j<<1|1] << k | code->x[j<<1|0]; + y = yak_hash64_64(y); + if (y % MODE_VALUE > 3) return 0; + int s = 64 - k; + uint64_t z = code->x[j<<1|1] >> s ^ y << (s + s) >> (s + s); + *get_sub_ID = y >> suffix_bits | z << (64 - suffix_bits); + *get_sub_key = y & suffix_mode; + return 1; +} + +inline int get_sub_table(uint64_t* get_sub_ID, uint64_t* get_sub_key, uint64_t suffix_mode, int suffix_bits, Hash_code* code, int k) +{ // not really working for k<=32 + return ha_get_sub_table_long(get_sub_ID, get_sub_key, suffix_mode, suffix_bits, code, k); +} inline int insert_Total_Count_Table(Total_Count_Table* TCB, Hash_code* code, int k) { @@ -286,7 +268,6 @@ inline int insert_Total_Count_Table(Total_Count_Table* TCB, Hash_code* code, int khint_t t; int absent; - while (__sync_lock_test_and_set(&TCB->sub_h_lock[sub_ID].lock, 1)) { while (TCB->sub_h_lock[sub_ID].lock); @@ -310,7 +291,6 @@ inline int insert_Total_Count_Table(Total_Count_Table* TCB, Hash_code* code, int inline int get_Total_Count_Table(Total_Count_Table* TCB, Hash_code* code, int k) { - uint64_t sub_ID, sub_key; if(!get_sub_table(&sub_ID, &sub_key, TCB->suffix_mode, TCB->suffix_bits, code, k)) { @@ -319,7 +299,6 @@ inline int get_Total_Count_Table(Total_Count_Table* TCB, Hash_code* code, int k) khint_t t; - ///query hash table，key is k t = ha_ct_get(TCB->sub_h[sub_ID], sub_key); if (t != kh_end(TCB->sub_h[sub_ID])) @@ -330,15 +309,10 @@ inline int get_Total_Count_Table(Total_Count_Table* TCB, Hash_code* code, int k) { return 0; } - - } - - inline uint64_t get_Total_Pos_Table(Total_Pos_Table* PCB, Hash_code* code, int k, uint64_t* r_sub_ID) { - uint64_t sub_ID, sub_key; if(!get_sub_table(&sub_ID, &sub_key, PCB->suffix_mode, PCB->suffix_bits, code, k)) { @@ -347,7 +321,6 @@ inline uint64_t get_Total_Pos_Table(Total_Pos_Table* PCB, Hash_code* code, int k khint_t t; - ///query hash table，key is k t = ha_pt_get(PCB->sub_h[sub_ID], sub_key); if (t != kh_end(PCB->sub_h[sub_ID])) @@ -359,8 +332,6 @@ inline uint64_t get_Total_Pos_Table(Total_Pos_Table* PCB, Hash_code* code, int k { return (uint64_t)-1; } - - } inline uint64_t count_Total_Pos_Table(Total_Pos_Table* PCB, Hash_code* code, int k) @@ -405,7 +376,6 @@ inline uint64_t insert_Total_Pos_Table(Total_Pos_Table* PCB, Hash_code* code, in if (occ) { - while (__sync_lock_test_and_set(&PCB->sub_h_lock[sub_ID].lock, 1)) { while (PCB->sub_h_lock[sub_ID].lock); @@ -413,16 +383,16 @@ inline uint64_t insert_Total_Pos_Table(Total_Pos_Table* PCB, Hash_code* code, in if (list[0].offset + 1 < occ) { - list[0].offset++; + list[0].offset++; // if not the last k-mer, this field is reused to keep the number of inserted positions list[list[0].offset].readID = readID; - ///list[list[0].offset].readID = readID|direction; list[list[0].offset].offset = pos; + list[list[0].offset].rev = ha_code2rev(code); } - else + else // now comes to the last k-mer position; then save it to list[0] { list[0].readID = readID; - ///list[0].readID = readID|direction; list[0].offset = pos; + list[0].rev = ha_code2rev(code); flag = 1; } @@ -434,18 +404,14 @@ inline uint64_t insert_Total_Pos_Table(Total_Pos_Table* PCB, Hash_code* code, in qsort(list, occ, sizeof(k_mer_pos), cmp_k_mer_pos); } - return 1; } else { return 0; } - - } - void init_Total_Count_Table(int k, Total_Count_Table* TCB); void init_Total_Pos_Table(Total_Pos_Table* TCB, Total_Count_Table* pre_TCB); void destory_Total_Count_Table(Total_Count_Table* TCB); @@ -473,12 +439,6 @@ uint64_t n_end_pos, uint8_t n_direction); -void merge_k_mer_pos_list_alloc_heap_sort(k_mer_pos_list_alloc* list, Candidates_list* candidates, HeapSq* HBT); - -void Init_Heap(HeapSq* HBT); -void destory_Heap(HeapSq* HBT); -void clear_Heap(HeapSq* HBT); - void init_overlap_region_alloc(overlap_region_alloc* list); void clear_overlap_region_alloc(overlap_region_alloc* list); void destory_overlap_region_alloc(overlap_region_alloc* list); diff --git a/Makefile b/Makefile index d8e1773..8323c77 100644 --- a/Makefile +++ b/Makefile @@ -50,7 +50,7 @@ POA.o: POA.h Hash_Table.h khashl.h kmer.h Process_Read.h kseq.h Overlaps.h POA.o: kvec.h kdq.h CommandLines.h Correct.h Levenshtein_distance.h Process_Read.o: Process_Read.h kseq.h Overlaps.h kvec.h kdq.h CommandLines.h Trio.o: khashl.h kthread.h Process_Read.h kseq.h Overlaps.h kvec.h kdq.h -Trio.o: CommandLines.h Trio.h +Trio.o: CommandLines.h Trio.h kmer.h kmer.o: kmer.h Process_Read.h kseq.h Overlaps.h kvec.h kdq.h CommandLines.h kthread.o: kthread.h main.o: CommandLines.h Process_Read.h kseq.h Overlaps.h kvec.h kdq.h diff --git a/Trio.cpp b/Trio.cpp index d9d1657..4bf4b8b 100644 --- a/Trio.cpp +++ b/Trio.cpp @@ -8,6 +8,7 @@ #include "Process_Read.h" #include "Trio.h" #include "CommandLines.h" +#include "kmer.h" #define CALLOC(ptr, len) ((ptr) = (__typeof__(ptr))calloc((len), sizeof(*(ptr)))) #define MALLOC(ptr, len) ((ptr) = (__typeof__(ptr))malloc((len) * sizeof(*(ptr)))) @@ -37,30 +38,6 @@ unsigned char seq_nt4_table[256] = { // translate ACGT to 0123 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4 }; -static inline uint64_t yak_hash64(uint64_t key, uint64_t mask) // invertible integer hash function -{ - key = (~key + (key << 21)) & mask; // key = (key << 21) - key - 1; - key = key ^ key >> 24; - key = ((key + (key << 3)) + (key << 8)) & mask; // key * 265 - key = key ^ key >> 14; - key = ((key + (key << 2)) + (key << 4)) & mask; // key * 21 - key = key ^ key >> 28; - key = (key + (key << 31)) & mask; - return key; -} - -static inline uint64_t yak_hash64_64(uint64_t key) -{ - key = ~key + (key << 21); - key = key ^ key >> 24; - key = (key + (key << 3)) + (key << 8); - key = key ^ key >> 14; - key = (key + (key << 2)) + (key << 4); - key = key ^ key >> 28; - key = key + (key << 31); - return key; -} - static inline uint64_t yak_hash_long(uint64_t x[4]) { int j = x[1] < x[3]? 0 : 1; @@ -361,4 +338,4 @@ void trio_partition() fprintf(stderr, "Trio binning has been done.\n"); fprintf(stderr, "%-30s%18.2f\n\n", "Trio binning time:", Get_T() - start_time); -} \ No newline at end of file +} diff --git a/kmer.cpp b/kmer.cpp index 0d234cf..811f8ec 100644 --- a/kmer.cpp +++ b/kmer.cpp @@ -10,14 +10,6 @@ void init_HPC_seq(HPC_seq* seq, char* str, long long l) seq->str = str; } -void init_Hash_code(Hash_code* code) -{ - code->x[0] = 0; - code->x[1] = 0; -} - - - void init_small_hash_table(small_hash_table* x) { x->size = 0; diff --git a/kmer.h b/kmer.h index 016a531..6462f57 100644 --- a/kmer.h +++ b/kmer.h @@ -16,7 +16,7 @@ typedef struct { //can represent at most 64-mer - uint64_t x[2]; + uint64_t x[4]; } Hash_code; typedef struct { @@ -83,16 +83,45 @@ inline uint64_t get_HPC_code(HPC_seq* seq, uint64_t* end_pos) } -inline void k_mer_append(Hash_code* code, uint64_t c, int k) +inline void init_Hash_code(Hash_code* code) { - - uint64_t mask = ALL >> (64 -k); - - code->x[0] = ((code->x[0]<<1) | (c&1)) & mask; - code->x[1] = ((code->x[1]<<1) | (c>>1)) & mask; + code->x[0] = code->x[1] = code->x[2] = code->x[3] = 0; } -inline void Hashcode_to_string(Hash_code* code, char* str, int k) +inline void k_mer_append(Hash_code* code, uint64_t c, int k) +{ + uint64_t mask = ALL >> (64 - k), shift = k - 1; + code->x[0] = ((code->x[0]<<1) | (c&1)) & mask; + code->x[1] = ((code->x[1]<<1) | (c>>1)) & mask; + code->x[2] = code->x[2] >> 1 | (uint64_t)(1 - (c&1)) << shift; + code->x[3] = code->x[3] >> 1 | (uint64_t)(1 - (c>>1)) << shift; +} + +static inline uint64_t yak_hash64(uint64_t key, uint64_t mask) // invertible integer hash function +{ + key = (~key + (key << 21)) & mask; // key = (key << 21) - key - 1; + key = key ^ key >> 24; + key = ((key + (key << 3)) + (key << 8)) & mask; // key * 265 + key = key ^ key >> 14; + key = ((key + (key << 2)) + (key << 4)) & mask; // key * 21 + key = key ^ key >> 28; + key = (key + (key << 31)) & mask; + return key; +} + +static inline uint64_t yak_hash64_64(uint64_t key) +{ + key = ~key + (key << 21); + key = key ^ key >> 24; + key = (key + (key << 3)) + (key << 8); + key = key ^ key >> 14; + key = (key + (key << 2)) + (key << 4); + key = key ^ key >> 28; + key = key + (key << 31); + return key; +} + +inline void Hashcode_to_string(Hash_code* code, char* str, int k) // FIXME: not working { uint8_t c; int i; @@ -108,7 +137,5 @@ inline void Hashcode_to_string(Hash_code* code, char* str, int k) } void init_HPC_seq(HPC_seq* seq, char* str, long long l); -void init_Hash_code(Hash_code* code); - -#endif \ No newline at end of file +#endif