diff --git a/.vscode/settings.json b/.vscode/settings.json
index 44c3105..a4dc91c 100644
--- a/.vscode/settings.json
+++ b/.vscode/settings.json
@@ -10,6 +10,11 @@
         "*.tcc": "cpp",
         "bitset": "cpp",
         "algorithm": "cpp",
-        "hashtable": "cpp"
+        "hashtable": "cpp",
+        "string_view": "cpp",
+        "list": "cpp",
+        "string": "cpp",
+        "array": "cpp",
+        "utility": "cpp"
     }
 }
\ No newline at end of file
diff --git a/Assembly.cpp b/Assembly.cpp
index 2429e52..fea217d 100644
--- a/Assembly.cpp
+++ b/Assembly.cpp
@@ -21,6 +21,7 @@ long long total_potiental_matched_overlap_0 = 0;
 long long total_potiental_matched_overlap_1 = 0;
 long long total_num_read_base = 0;
 long long total_num_correct_base = 0;
+int roundID = 0;
 
 
 long long complete_threads = 0;
@@ -139,12 +140,146 @@ void* Perform_Counting(void* arg)
     }
 
     destory_R_buffer_block(&curr_sub_block);
-    free(arg);
+    ///free(arg);
 
 }
 
 
 
+void* Perform_Counting_non_first(void* arg)
+{
+    int thr_ID = *((int*)arg);
+
+    int i = 0;
+    HPC_seq HPC_read;
+
+
+    long long read_number = 0;
+    long long select_k_mer_number = 0 ;
+    long long k_mer_number = 0 ;
+
+	int file_flag = 1;
+
+    uint64_t code;
+
+    uint64_t end_pos;
+
+    Hash_code k_code;
+
+    int avalible_k = 0;
+
+    UC_Read g_read;
+    init_UC_Read(&g_read);
+
+    for (i = thr_ID; i < R_INF.total_reads; i = i + thread_num)
+    {
+        recover_UC_Read(&g_read, &R_INF, i);
+        ///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,k_mer_length);
+                avalible_k++;
+                if (avalible_k>=k_mer_length)
+                {
+                    ///插入
+                    if(insert_Total_Count_Table(&TCB, &k_code, k_mer_length))
+                    {
+                        select_k_mer_number++;
+                    }
+
+                    k_mer_number++;
+                    
+                }
+                
+            }
+            else
+            {
+                avalible_k = 0;
+                init_Hash_code(&k_code);
+            }
+            
+        }
+    }
+
+    destory_UC_Read(&g_read);
+}
+
+
+
+void* Build_hash_table_non_first(void* arg)
+{
+    int thr_ID = *((int*)arg);
+
+    int i = 0;
+    HPC_seq HPC_read;
+
+	int file_flag = 1;
+
+    uint64_t code;
+    uint64_t end_pos;
+
+    ///long long HPC_base;
+
+    Hash_code k_code;
+
+    int avalible_k = 0;
+
+    UC_Read g_read;
+    init_UC_Read(&g_read);
+
+    for (i = thr_ID; i < R_INF.total_reads; i = i + thread_num)
+    {
+        recover_UC_Read(&g_read, &R_INF, i);
+        ///forward strand
+        init_HPC_seq(&HPC_read, g_read.seq, g_read.length);
+        init_Hash_code(&k_code);
+
+        avalible_k = 0;
+
+        ///HPC_base = 0;
+
+        while ((code = get_HPC_code(&HPC_read, &end_pos)) != 6)
+        {
+            if(code < 4)
+            {
+                k_mer_append(&k_code,code,k_mer_length);
+                avalible_k++;
+                if (avalible_k>=k_mer_length)
+                {
+
+                    ///选取的k-mer满足两个要求
+                    ///1. hash(k-mer) % 101 <= 3
+                    ///2. occ(k-mer)要满足范围
+                    ///TCB表中的元素仅满足第一个要求，而PCB表中的元素满足两个要求
+                    ///所以如果当前k-mer在PCB表中存在，则他的位置一定要加入到候选位置中去
+                    ///insert_Total_Pos_Table(&PCB, &k_code, k_mer_length, curr_sub_block.read[i].ID, HPC_base - k_mer_length + 1, FORWARD);
+                    insert_Total_Pos_Table(&PCB, &k_code, k_mer_length, 
+                        i, end_pos);
+                }
+                
+            }
+            else
+            {
+                avalible_k = 0;
+                init_Hash_code(&k_code);
+            }
+
+            ///HPC_base++;
+        }
+    }
+
+    destory_UC_Read(&g_read);
+}
+
+
+
 
 
 
@@ -275,7 +410,7 @@ void* Build_hash_table(void* arg)
     }
 
     destory_R_buffer_block(&curr_sub_block);
-    free(arg);
+    ///free(arg);
 
 
 
@@ -296,20 +431,21 @@ void Counting_multiple_thr()
 
     fprintf(stdout, "Begin Counting ...... \n");
 
-    init_kseq(read_file_name);
-
-    init_All_reads(&R_INF);
-
     init_Total_Count_Table(k_mer_length, &TCB);
 
     pthread_t inputReadsHandle;
-    
-    init_R_buffer(thread_num);
 
     int *is_insert = (int*)malloc(sizeof(*is_insert));
 	*is_insert = 1;
 
-    pthread_create(&inputReadsHandle, NULL, input_reads_muti_threads, (void*)is_insert);
+    if (roundID == 0)
+    {
+        init_kseq(read_file_name);
+        init_All_reads(&R_INF);
+        init_R_buffer(thread_num);
+        pthread_create(&inputReadsHandle, NULL, input_reads_muti_threads, (void*)is_insert);
+    }
+    
 
     pthread_t *_r_threads;
 
@@ -321,14 +457,19 @@ void Counting_multiple_thr()
 	{
 		int *arg = (int*)malloc(sizeof(*arg));
 		*arg = i;
-
-		pthread_create(_r_threads + i, NULL, Perform_Counting, (void*)arg);
-
+        if (roundID == 0)
+        {
+            pthread_create(_r_threads + i, NULL, Perform_Counting, (void*)arg);
+        }
+        else
+        {
+            pthread_create(_r_threads + i, NULL, Perform_Counting_non_first, (void*)arg);
+        }
 	}
 
 
 
-    pthread_join(inputReadsHandle, NULL);
+    
 
     for (i = 0; i<thread_num; i++)
 		pthread_join(_r_threads[i], NULL);
@@ -338,8 +479,11 @@ void Counting_multiple_thr()
     ///destory_R_buffer();
     
     ///destory_Total_Count_Table(&TCB);
-
-    destory_kseq();
+    if (roundID == 0)
+    {
+        pthread_join(inputReadsHandle, NULL);
+        destory_kseq();
+    }
 
     fprintf(stdout, "Finish Counting ...... \n");
 
@@ -368,21 +512,24 @@ void Build_hash_table_multiple_thr()
     fprintf(stdout, "%-30s%18.2f\n\n", "Traverse time:", Get_T() - T_start_time);
     fflush(stdout);
    
-
-
-    init_kseq(read_file_name);
-
-    clear_R_buffer();
-
-    malloc_All_reads(&R_INF);
+    ///at this moment, TCB can be free
+    destory_Total_Count_Table(&TCB);
 
     pthread_t inputReadsHandle;
 
     int *is_insert = (int*)malloc(sizeof(*is_insert));
 	*is_insert = 0;
 
-    pthread_create(&inputReadsHandle, NULL, input_reads_muti_threads, (void*)is_insert);
+    if (roundID == 0)
+    {
+        init_kseq(read_file_name);
 
+        clear_R_buffer();
+
+        malloc_All_reads(&R_INF);
+
+        pthread_create(&inputReadsHandle, NULL, input_reads_muti_threads, (void*)is_insert);
+    }
 
 
 
@@ -397,36 +544,49 @@ void Build_hash_table_multiple_thr()
 		int *arg = (int*)malloc(sizeof(*arg));
 		*arg = i;
 
-		pthread_create(_r_threads + i, NULL, Build_hash_table, (void*)arg);
-
+        if (roundID == 0)
+        {
+            pthread_create(_r_threads + i, NULL, Build_hash_table, (void*)arg);
+        }
+        else
+        {
+            pthread_create(_r_threads + i, NULL, Build_hash_table_non_first, (void*)arg);
+        }
 	}
 
-    pthread_join(inputReadsHandle, NULL);
+    if (roundID == 0)
+    {
+        pthread_join(inputReadsHandle, NULL);
+    }
 
     for (i = 0; i<thread_num; i++)
 		pthread_join(_r_threads[i], NULL);
 
     free(_r_threads);
 
-    destory_kseq();
+    
 
     fprintf(stdout, "Finish Building hash table ...... \n");
 
     fprintf(stdout, "%-30s%18.2f\n\n", "Build hash table time:", Get_T() - start_time);
 
-    ///destory_Total_Count_Table(&TCB);
-    if (write_index_to_disk)
+    if (roundID == 0)
     {
-        write_Total_Pos_Table(&PCB, read_file_name);
-        ///destory_Total_Pos_Table(&PCB);
-        ///load_Total_Pos_Table(&PCB, read_file_name);
-
-        write_All_reads(&R_INF, read_file_name);
-        ///destory_All_reads(&R_INF);
-        ///load_All_reads(&R_INF, read_file_name);
+        destory_kseq();
+        destory_R_buffer();
+        ///destory_Total_Count_Table(&TCB);
+        if (write_index_to_disk)
+        {
+            write_Total_Pos_Table(&PCB, read_file_name);
+            ///destory_Total_Pos_Table(&PCB);
+            ///load_Total_Pos_Table(&PCB, read_file_name);
+            write_All_reads(&R_INF, read_file_name);
+            ///destory_All_reads(&R_INF);
+            ///load_All_reads(&R_INF, read_file_name);
+        }
     }
 
-    destory_Total_Count_Table(&TCB);
+    ///destory_Total_Count_Table(&TCB);
     
     free(is_insert);
     
@@ -787,6 +947,294 @@ int get_required_read(const char *required_name, long long RID, All_reads* R_INF
 }
 
 
+void get_corrected_read_from_cigar(Cigar_record* cigar, char* pre_read, int pre_length,
+char* new_read, int* new_length)
+{
+    int i, j;
+    int pre_i, new_i;
+    int operation, operation_length;
+    pre_i = new_i = 0;
+    int diff_char_i = 0;
+
+
+    for (i = 0; i < cigar->length; i++)
+    {
+        operation = Get_Cigar_Type(cigar->record[i]);
+        operation_length = Get_Cigar_Length(cigar->record[i]);
+
+        if (operation == 0)
+        {
+            memcpy(new_read + new_i, pre_read + pre_i, operation_length);
+            pre_i = pre_i + operation_length;
+            new_i = new_i + operation_length;
+        }
+        else if (operation == 1)
+        {
+            
+            for (j = 0; j < operation_length; j++)
+            {
+                new_read[new_i] = Get_MisMatch_Base(cigar->lost_base[diff_char_i]);
+                new_i++;
+                diff_char_i++;
+            }
+            pre_i = pre_i + operation_length;
+            
+            /**
+            pre_i = pre_i + operation_length;
+            memcpy(new_read + new_i, cigar->lost_base + diff_char_i, operation_length);
+            new_i = new_i + operation_length;
+            diff_char_i = diff_char_i + operation_length;
+            **/
+            
+        }
+        else if (operation == 3)
+        {
+            pre_i = pre_i + operation_length;
+            diff_char_i = diff_char_i + operation_length;
+        }
+        else if (operation == 2)
+        {
+            memcpy(new_read + new_i, cigar->lost_base + diff_char_i, operation_length);
+            new_i = new_i + operation_length;
+            diff_char_i = diff_char_i + operation_length;
+        }
+    }
+    *new_length = new_i;
+}
+
+
+void get_uncorrected_read_from_cigar(Cigar_record* cigar, char* new_read, int new_length, char* pre_read, int* pre_length)
+{
+    int i, j;
+    int pre_i, new_i;
+    int operation, operation_length;
+    pre_i = new_i = 0;
+    int diff_char_i = 0;
+
+
+    for (i = 0; i < cigar->length; i++)
+    {
+        operation = Get_Cigar_Type(cigar->record[i]);
+        operation_length = Get_Cigar_Length(cigar->record[i]);
+
+
+        if (operation == 0)
+        {
+            memcpy(pre_read + pre_i, new_read + new_i, operation_length);
+            pre_i = pre_i + operation_length;
+            new_i = new_i + operation_length;
+        }
+        else if (operation == 1)
+        {
+            
+            for (j = 0; j < operation_length; j++)
+            {
+                pre_read[pre_i] = Get_Match_Base(cigar->lost_base[diff_char_i]);
+                pre_i++;
+                diff_char_i++;
+            }
+            new_i = new_i + operation_length;
+        }
+        else if (operation == 3)
+        {
+            memcpy(pre_read + pre_i, cigar->lost_base + diff_char_i, operation_length);
+            pre_i = pre_i + operation_length;
+            diff_char_i = diff_char_i + operation_length;
+        }
+        else if (operation == 2)
+        {
+            new_i = new_i + operation_length;
+            diff_char_i = diff_char_i + operation_length;
+        }
+    }
+
+    *pre_length = pre_i;
+}
+
+int debug_cigar(Cigar_record* cigar, char* pre_read, int pre_length,
+char* new_read, int new_length, int correct_base)
+{
+    int i;
+    int total_errors = 0;
+    for (i = 0; i < cigar->length; i++)
+    {
+        if (Get_Cigar_Type(cigar->record[i]) > 0)
+        {
+            total_errors = total_errors + Get_Cigar_Length(cigar->record[i]);
+        }
+    }
+
+    if(total_errors!=correct_base)
+    {
+        fprintf(stderr, "total_errors: %d, correct_base: %d\n", total_errors, correct_base);
+    }
+
+    int pre_i, new_i;
+    int operation, operation_length;
+    pre_i = new_i = 0;
+
+
+    for (i = 0; i < cigar->length; i++)
+    {
+        operation = Get_Cigar_Type(cigar->record[i]);
+        operation_length = Get_Cigar_Length(cigar->record[i]);
+
+        if (operation == 0)
+        {
+            pre_i = pre_i + operation_length;
+            new_i = new_i + operation_length;
+        }
+
+        if (operation == 1)
+        {
+            pre_i = pre_i + operation_length;
+            new_i = new_i + operation_length;
+        }
+
+        if (operation == 3)
+        {
+            pre_i = pre_i + operation_length;
+        }
+
+        if (operation == 2)
+        {
+            new_i = new_i + operation_length;
+        }
+        
+    }
+
+    /**
+    fprintf(stderr, "total_errors: %d, correct_base: %d, length: %d, lost_base_length: %d\n", 
+    total_errors, correct_base, cigar->length, cigar->lost_base_length);
+    **/
+
+    if (pre_i != pre_length || new_i != new_length)
+    {
+        fprintf(stderr, "pre_i: %d, pre_length: %d\n", pre_i, pre_length);
+        fprintf(stderr, "new_i: %d, new_length: %d\n", new_i, new_length);
+    }
+
+    char* tmp_seq = (char*)malloc(new_length + pre_length);
+    int tmp_length;
+
+    get_corrected_read_from_cigar(cigar, pre_read, pre_length, tmp_seq, &tmp_length);
+
+    if(tmp_length != new_length)
+    {
+        fprintf(stderr, "tmp_length: %d, new_length: %d\n", tmp_length, new_length);
+    }
+
+    if(memcmp(new_read, tmp_seq, new_length)!=0)
+    {
+        fprintf(stderr, "error new string\n");
+    }
+
+
+
+
+    
+    get_uncorrected_read_from_cigar(cigar, new_read, new_length, tmp_seq, &tmp_length);
+
+    
+    if(tmp_length != pre_length)
+    {
+        fprintf(stderr, "tmp_length: %d, pre_length: %d\n", tmp_length, pre_length);
+    }
+
+    if(memcmp(pre_read, tmp_seq, pre_length)!=0)
+    {
+        fprintf(stderr, "error pre string\n");
+    }
+    
+
+    free(tmp_seq);
+
+    if(cigar->new_read_length != new_length)
+    {
+        fprintf(stderr, "cigar->new_read_length: %d, new_length: %d\n", cigar->new_read_length, new_length);
+    }
+    
+    
+    
+}
+
+inline void push_cigar(Compressed_Cigar_record* records, long long ID, Cigar_record* input)
+{
+
+    if (input->length > records[ID].size)
+    {
+        records[ID].size = input->length;
+        records[ID].record = (uint32_t*)realloc(records[ID].record, records[ID].size*sizeof(uint32_t));     
+    }
+    records[ID].length = input->length;
+    memcpy(records[ID].record, input->record, input->length*sizeof(uint32_t));
+
+    if (input->lost_base_length > records[ID].lost_base_size)
+    {
+        records[ID].lost_base_size = input->lost_base_length;
+        records[ID].lost_base = (char*)realloc(records[ID].lost_base, records[ID].lost_base_size);
+    }
+    records[ID].lost_base_length = input->lost_base_length;
+    memcpy(records[ID].lost_base, input->lost_base, input->lost_base_length);
+
+    records[ID].new_length = input->new_read_length;
+    
+    
+}
+
+void just_debug(overlap_region_alloc* overlap_list, All_reads* R_INF)
+{
+    int j;
+    int i;
+    int y_id;
+    int y_strand;
+    int y_readLen;
+    int overlap_length;
+    int high_quality_overlaps = 0;
+    UC_Read g_read;
+    init_UC_Read(&g_read);
+
+    for (j = 0; j < overlap_list->length; j++)
+    {
+        y_id = overlap_list->list[j].y_id;
+        y_strand = overlap_list->list[j].y_pos_strand;
+        y_readLen = Get_READ_LENGTH((*R_INF), y_id);
+        overlap_length = overlap_list->list[j].x_pos_e - overlap_list->list[j].x_pos_s + 1;
+        if (overlap_length * OVERLAP_THRESHOLD <=  overlap_list->list[j].align_length)
+        {
+            high_quality_overlaps++;
+
+            fprintf(stderr, "i: %d, overlap_length: %d, align_length: %d, y_strand: %d\n", 
+                    high_quality_overlaps, overlap_length, overlap_list->list[j].align_length, y_strand);
+            fprintf(stderr, "x_pos_s: %d, x_pos_e: %d\n", 
+                    overlap_list->list[j].x_pos_s, overlap_list->list[j].x_pos_e);
+
+            fprintf(stderr, "%.*s\n\n", Get_NAME_LENGTH((*R_INF), y_id), Get_NAME((*R_INF), y_id));
+
+        }
+    }
+
+
+    recover_UC_Read_RC(&g_read, R_INF, overlap_list->list[0].x_id);
+
+    for (i = 0; i < g_read.length && i < 81; i++)
+    {
+        fprintf(stderr, "%c", g_read.seq[i]);
+    }
+
+    fprintf(stderr, "\n");
+
+    fprintf(stderr, "x_length: %d\n", g_read.length);
+
+
+    
+
+    fprintf(stderr, "high_quality_overlaps: %d, overlap_list->length: %d\n", high_quality_overlaps,overlap_list->length);
+
+    destory_UC_Read(&g_read);
+
+}
+
 void* Overlap_calculate_heap_merge(void* arg)
 {
     /************需要注释掉**********/
@@ -850,17 +1298,19 @@ void* Overlap_calculate_heap_merge(void* arg)
 
     init_buffer_sub_block(&current_sub_buffer);
 
-    
+    Cigar_record current_cigar;
+    init_Cigar_record(&current_cigar);
+
+    haplotype_evdience_alloc hap;
+    InitHaplotypeEvdience(&hap);
 
     for (i = thr_ID; i < R_INF.total_reads; i = i + thread_num)
     {
-        /**
-        if(!get_required_read("m54238_180916_191625\/43253934\/ccs", i, &R_INF))
-        {
-            continue;
-        }
-        **/
         
+        
+        
+
+        clear_Cigar_record(&current_cigar);
         clear_Heap(&heap);
         clear_Candidates_list(&l);
         ///clear_Candidates_list(&debug_l);
@@ -996,14 +1446,38 @@ void* Overlap_calculate_heap_merge(void* arg)
        ///clear_Graph(&POA_Graph);
 
        correct_overlap(&overlap_list, &R_INF, &g_read, &correct, &overlap_read, &POA_Graph,
-       &matched_overlap_0, &matched_overlap_1, &potiental_matched_overlap_0, &potiental_matched_overlap_1);
+       &matched_overlap_0, &matched_overlap_1, &potiental_matched_overlap_0, &potiental_matched_overlap_1,
+       &current_cigar, &hap);
       
         num_read_base = num_read_base + g_read.length;
         num_correct_base = num_correct_base + correct.corrected_base;
 
-        output_read_to_buffer(i, &R_INF, correct.corrected_read, correct.corrected_read_length, &current_sub_buffer);
-        ///output_read_to_buffer(i, &R_INF, g_read.seq, g_read.length, &current_sub_buffer);
+        push_cigar(R_INF.cigars, i, &current_cigar);
 
+        /**
+        if(memcmp("m54334_180926_225337/39780640/ccs", Get_NAME(R_INF, i), Get_NAME_LENGTH(R_INF, i)) == 0)
+        {
+            fprintf(stderr, "i: %d\n", i);
+            ///fprintf(stderr, "overlap_list.length: %d\n", overlap_list.length);
+            just_debug(&overlap_list, &R_INF);
+        }
+        **/
+        
+        
+
+
+
+        ///output_read_to_buffer(i, &R_INF, correct.corrected_read, correct.corrected_read_length, &current_sub_buffer);
+        /**
+        debug_cigar(&current_cigar, g_read.seq, g_read.length, correct.corrected_read, correct.corrected_read_length, 
+        correct.corrected_base);
+        **/
+        
+        
+        
+
+
+        
         
         /**
         POA_i = 0;
@@ -1109,10 +1583,12 @@ void* Overlap_calculate_heap_merge(void* arg)
     destory_Graph(&POA_Graph);
     destory_UC_Read(&g_read);
     destory_UC_Read(&overlap_read);
-    
+    destory_Cigar_record(&current_cigar);
 
     destory_Correct_dumy(&correct);
 
+    destoryHaplotypeEvdience(&hap);
+
 
     pthread_mutex_lock(&statistics);
     total_matched_overlap_0 += matched_overlap_0;
@@ -1133,9 +1609,115 @@ void* Overlap_calculate_heap_merge(void* arg)
         fprintf(stderr, "total_num_correct_base: %llu\n", total_num_correct_base);
     }
 	pthread_mutex_unlock(&statistics);
+
+    free(arg);
 }
 
 
+void* Save_corrected_reads(void* arg)
+{
+    int thr_ID = *((int*)arg);
+    long long i, j;
+    UC_Read g_read;
+    init_UC_Read(&g_read);
+
+    int new_read_size = 10000;
+    char* new_read = (char*)malloc(new_read_size);
+    Cigar_record cigar;
+    int new_read_length;
+    uint64_t N_occ;
+    
+    for (i = thr_ID; i < R_INF.total_reads; i = i + thread_num)
+    {
+        recover_UC_Read(&g_read, &R_INF, i);
+
+        if(R_INF.cigars[i].new_length>new_read_size)
+        {
+            new_read_size = R_INF.cigars[i].new_length;
+            new_read = (char*)realloc(new_read, new_read_size);
+        }
+
+        cigar.length = R_INF.cigars[i].length;
+        cigar.lost_base_length = R_INF.cigars[i].lost_base_length;
+        cigar.record = R_INF.cigars[i].record;
+        cigar.lost_base = R_INF.cigars[i].lost_base;
+
+        get_corrected_read_from_cigar(&cigar, g_read.seq, g_read.length, new_read, &new_read_length);
+
+        ///need modification
+        reverse_complement(new_read, new_read_length);
+
+
+        N_occ = 0;
+        for (j = 0; j < new_read_length; j++)
+        {
+            if(new_read[j] == 'N')
+            {
+                N_occ++;
+            }
+        }
+        
+        if(R_INF.read_size[i] < new_read_length)
+        {
+            R_INF.read_size[i] = new_read_length;
+            R_INF.read_sperate[i] = (uint8_t*)realloc(R_INF.read_sperate[i], R_INF.read_size[i]/4+1);
+        }
+
+        R_INF.read_length[i] = new_read_length;
+        
+
+        compress_base(Get_READ(R_INF, i),
+            new_read, new_read_length, 
+            &R_INF.N_site[i], N_occ);
+        
+        
+    }
+
+    destory_UC_Read(&g_read);
+    free(new_read);
+    free(arg);
+}
+
+
+void Output_corrected_reads()
+{
+    long long i, j;
+    UC_Read g_read;
+    init_UC_Read(&g_read);
+    FILE* output_file = fopen(output_file_name, "w");
+
+
+    if(number_of_round % 2 == 0)
+    {
+        for (i = 0; i < R_INF.total_reads; i++)
+        {
+            recover_UC_Read(&g_read, &R_INF, i);
+            fwrite(">", 1, 1, output_file);
+            fwrite(Get_NAME(R_INF, i), 1, Get_NAME_LENGTH(R_INF, i), output_file);
+            fwrite("\n", 1, 1, output_file);
+            fwrite(g_read.seq, 1, g_read.length, output_file);
+            fwrite("\n", 1, 1, output_file);
+        }
+    }
+    else
+    {
+        for (i = 0; i < R_INF.total_reads; i++)
+        {
+            recover_UC_Read_RC(&g_read, &R_INF, i);
+            fwrite(">", 1, 1, output_file);
+            fwrite(Get_NAME(R_INF, i), 1, Get_NAME_LENGTH(R_INF, i), output_file);
+            fwrite("\n", 1, 1, output_file);
+            fwrite(g_read.seq, 1, g_read.length, output_file);
+            fwrite("\n", 1, 1, output_file);
+        }
+    }
+    
+    
+
+
+    destory_UC_Read(&g_read);
+    fclose(output_file);
+}
 
 
 
@@ -1145,20 +1727,19 @@ void* Overlap_calculate_heap_merge(void* arg)
 void Overlap_calculate_multipe_thr()
 {
 
-    fprintf(stderr, "k_mer_min_freq: %lld, CORRECT_THRESHOLD: %f\n", k_mer_min_freq, CORRECT_THRESHOLD);
+
     
     double start_time = Get_T();
-
-    init_output_buffer(thread_num);
-
     pthread_t outputResultSinkHandle;
+    /**
+    if (roundID == number_of_round - 1)
+    {
+        init_output_buffer(thread_num);
+        pthread_create(&outputResultSinkHandle, NULL, pop_buffer, NULL);
+    }
+    **/
 
-    pthread_create(&outputResultSinkHandle, NULL, pop_buffer, NULL);
-
-
-
-    fprintf(stdout, "R_INF.total_reads: %llu\n", R_INF.total_reads);
-
+    
 
     fprintf(stdout, "Begin Overlap Calculate ...... \n");
 
@@ -1174,7 +1755,6 @@ void Overlap_calculate_multipe_thr()
 		*arg = i;
 
 		pthread_create(_r_threads + i, NULL, Overlap_calculate_heap_merge, (void*)arg);
-        //pthread_create(_r_threads + i, NULL, Overlap_calculate, (void*)arg);
 
 	}
     
@@ -1182,17 +1762,55 @@ void Overlap_calculate_multipe_thr()
     for (i = 0; i<thread_num; i++)
 		pthread_join(_r_threads[i], NULL);
 
-    pthread_join(outputResultSinkHandle, NULL);
+    /**
+    if (roundID == number_of_round - 1)
+    {
+        pthread_join(outputResultSinkHandle, NULL);
+        destory_output_buffer();
+        ///destory_All_reads(&R_INF);
+    }
+    **/
+    
 
     free(_r_threads);
 
-    destory_output_buffer();
 
-    
+    destory_Total_Pos_Table(&PCB);
 
     fprintf(stdout, "Finish Overlap Calculate.\n");
 
     fprintf(stdout, "%-30s%18.2f\n\n", "Calculate Overlap time:", Get_T() - start_time);    
+
+
+    start_time = Get_T();
+    _r_threads = (pthread_t *)malloc(sizeof(pthread_t)*thread_num);
+    for (i = 0; i < thread_num; i++)
+	{
+		int *arg = (int*)malloc(sizeof(*arg));
+		*arg = i;
+
+		pthread_create(_r_threads + i, NULL, Save_corrected_reads, (void*)arg);
+        //pthread_create(_r_threads + i, NULL, Overlap_calculate, (void*)arg);
+
+	}
+    for (i = 0; i<thread_num; i++)
+		pthread_join(_r_threads[i], NULL);
+    free(_r_threads);
+    fprintf(stdout, "%-30s%18.2f\n\n", "Save corrected read time:", Get_T() - start_time);   
+    ///fflush(stdout);
+
+    
+    ///only the last round can output read to disk
+    if (roundID == number_of_round - 1)
+    {
+        start_time = Get_T();
+        
+        Output_corrected_reads();
+        
+        fprintf(stdout, "%-30s%18.2f\n\n", "Output time:", Get_T() - start_time);    
+    }
+    
+
 }
 
 
@@ -1438,6 +2056,50 @@ void verify_Position_hash_table()
 }
 
 
+void Correct_Reads(int last_round)
+{
+    complete_threads = 0;
+    total_matched_overlap_0 = 0;
+    total_matched_overlap_1 = 0;
+    total_potiental_matched_overlap_0 = 0;
+    total_potiental_matched_overlap_1 = 0;
+    total_num_read_base = 0;
+    total_num_correct_base = 0;
+
+
+    if(last_round == 0)
+        return;
+
+    roundID = number_of_round - last_round;
+    fprintf(stdout, "Error correction: start the %d-th round ...\n", roundID);
+    ///only the first round correction can load index from disk
+    if (roundID == 0 & load_index_from_disk && load_pre_cauculated_index())
+    {
+        ;
+    }
+    else
+    {
+        Counting_multiple_thr();
+
+        Build_hash_table_multiple_thr();
+    }
+
+    
+    fprintf(stdout, "Total pos in hash tabe: %d\n", PCB.total_occ);
+    fprintf(stdout, "k_mer_min_freq in hashtable: %d\n", k_mer_min_freq);
+    fprintf(stdout, "k_mer_max_freq in hashtable: %d\n", k_mer_max_freq);
+
+    
+    
+    ///verify_Position_hash_table();
+    Overlap_calculate_multipe_thr();
+
+    fprintf(stdout, "Error correction: the %d-th round has been completed.\n", roundID);
+
+    Correct_Reads(last_round - 1);
+}
+
+
 
 
 
diff --git a/Assembly.h b/Assembly.h
index 0d853bc..6946f36 100644
--- a/Assembly.h
+++ b/Assembly.h
@@ -4,12 +4,15 @@
 #define FORWARD 0
 #define REVERSE_COMPLEMENT (0x8000000000000000)
 
-
+#define Get_Cigar_Type(RECORD) (RECORD&3)
+#define Get_Cigar_Length(RECORD) (RECORD>>2)
 
 void Counting_multiple_thr();
 void Build_hash_table_multiple_thr();
 int load_pre_cauculated_index();
 void Overlap_calculate_multipe_thr();
+void Correct_Reads(int last_round);
+
 /********************************for debug***************************************/
 void Verify_Counting();
 
diff --git a/CommandLines.cpp b/CommandLines.cpp
index 0a50fee..63c99a9 100644
--- a/CommandLines.cpp
+++ b/CommandLines.cpp
@@ -15,6 +15,7 @@ int k_mer_min_freq = 3;
 int k_mer_max_freq = 66;
 int load_index_from_disk = 0;
 int write_index_to_disk = 0;
+int number_of_round = 1;
 
 
 double Get_T(void)
@@ -42,19 +43,21 @@ int CommandLine_process (int argc, char *argv[])
         { "thread", ko_required_argument, 103},
         { "k_mer_min_freq", ko_required_argument, 104},
         { "k_mer_max_freq", ko_required_argument, 105},
+        { "round", ko_required_argument, 106},
         { NULL, 0, 0 }
     };
 
     ketopt_t opt = KETOPT_INIT;
 
     int i, c;
-    while ((c = ketopt(&opt, argc, argv, 1, "ht:o:q:k:lwm:n:", longopts)) >= 0) {
+    while ((c = ketopt(&opt, argc, argv, 1, "ht:o:q:k:lwm:n:r:", longopts)) >= 0) {
         if (c == 100 || c == 'h') Print_H();
         else if (c == 103 || c == 't') thread_num = atoi(opt.arg); 
         else if (c == 102 || c == 'o') output_file_name = opt.arg;
         else if (c == 101 || c == 'q') read_file_name = opt.arg;
         else if (c == 104 || c == 'n') k_mer_min_freq = atoi(opt.arg);
         else if (c == 105 || c == 'm') k_mer_max_freq = atoi(opt.arg);
+        else if (c == 106 || c == 'r') number_of_round = atoi(opt.arg);
         else if (c == 'k') k_mer_length = atoi(opt.arg); 
         else if (c == 'l') load_index_from_disk = 1; 
         else if (c == 'w') write_index_to_disk = 1; 
diff --git a/CommandLines.h b/CommandLines.h
index 0f85670..dd0525d 100644
--- a/CommandLines.h
+++ b/CommandLines.h
@@ -12,6 +12,7 @@ extern int k_mer_min_freq;
 extern int k_mer_max_freq;
 extern int load_index_from_disk;
 extern int write_index_to_disk;
+extern int number_of_round;
 
 
 int CommandLine_process (int argc, char *argv[]);
diff --git a/Correct.cpp b/Correct.cpp
index 09705ee..b37f147 100644
--- a/Correct.cpp
+++ b/Correct.cpp
@@ -445,6 +445,45 @@ long long o_len, int threashold, int error, long long* total_mis)
     edlibFreeAlignResult(result);
 }
 
+
+void fill_subregion(char* r, long long start_pos, long long length, uint8_t strand, All_reads* R_INF, long long ID, 
+int extra_begin, int extra_end)
+{
+    
+    recover_UC_Read_sub_region(r+extra_begin, start_pos, length, strand, R_INF, ID);
+    memset(r, 'N', extra_begin);
+    memset(r+extra_begin+length, 'N', extra_end);
+}
+
+void determine_overlap_region(int threshold, long long y_start, long long y_ID, long long Window_Len, All_reads* R_INF,
+int* r_extra_begin, int* r_extra_end, long long* r_y_start, long long* r_y_length)
+{
+    int extra_begin;
+    int extra_end;
+    long long currentIDLen;
+    long long o_len;
+
+    extra_begin = extra_end = 0;
+    ///y maybe less than 0
+    y_start = y_start - threshold;
+    ///the length of y
+    currentIDLen = Get_READ_LENGTH((*R_INF), y_ID);
+    o_len = MIN(Window_Len, currentIDLen - y_start);
+    extra_end = Window_Len - o_len;
+
+    if (y_start < 0)
+    {
+        extra_begin = -y_start;
+        y_start = 0;
+        o_len = o_len - extra_begin;
+    }
+
+    (*r_extra_begin) = extra_begin;
+    (*r_extra_end) = extra_end;
+    (*r_y_start) = y_start;
+    (*r_y_length) = o_len;
+}
+
 void verify_window(long long window_start, long long window_end, overlap_region_alloc* overlap_list,Correct_dumy* dumy, All_reads* R_INF,
 char* r_string)
 {
@@ -458,20 +497,20 @@ char* r_string)
     long long x_end, x_len;
     int end_site;
     unsigned int error;
-    /************需要注释掉********* */ 
-    // long long total_match=0;
-    // long long total_unmatch=0;
-    // long long total_mis=0;
-    /************需要注释掉********* */
     int groupLen = 0;
     int return_sites[GROUP_SIZE];
 	unsigned int return_sites_error[GROUP_SIZE];
     uint64_t overlapID[GROUP_SIZE];
     uint64_t y_startGroup[GROUP_SIZE];
+    int y_extra_begin[GROUP_SIZE];
+    int y_extra_end[GROUP_SIZE];
+    int extra_begin;
+    int extra_end;
 
     ///这些是整个window被完全覆盖的
     for (i = 0; i < dumy->length; i++)
     {
+        extra_begin = extra_end = 0;
         ///整个window被覆盖的话，read本身上的区间就是[window_start, window_end]
         x_len = WINDOW;
         currentID = dumy->overlapID[i];
@@ -479,173 +518,97 @@ char* r_string)
         ///y上的相对位置
         y_start = (x_start - overlap_list->list[currentID].x_pos_s) + overlap_list->list[currentID].y_pos_s;
 
-        // /************需要注释掉**********/
-        // if(x_start < overlap_list->list[currentID].x_pos_s)
-        // {
-        //     fprintf(stderr, "ERROR\n");
-        // }
-        // /************需要注释掉**********/
-        ///y上的起始
-        y_start = y_start - THRESHOLD;
-        if (y_start < 0)
-        {
-            y_start = 0;
-        }
-        ///当前y的长度
-        currentIDLen = Get_READ_LENGTH((*R_INF), overlap_list->list[currentID].y_id);
-        ///不能超过y的剩余长度
-        o_len = MIN(Window_Len, currentIDLen - y_start);
 
-        // /************需要注释掉**********/
-        // if(o_len < x_len)
-        // {
-        //     fprintf(stderr, "ERROR\n");
-        // }
-        // /************需要注释掉**********/
+  
+        determine_overlap_region(THRESHOLD, y_start, overlap_list->list[currentID].y_id, Window_Len, R_INF,
+        &extra_begin, &extra_end, &y_start, &o_len);
+
+        fill_subregion(dumy->overlap_region_group[groupLen], y_start, o_len, overlap_list->list[currentID].y_pos_strand, 
+        R_INF, overlap_list->list[currentID].y_id, extra_begin, extra_end);
+
+        y_extra_begin[groupLen] = extra_begin;
+        y_extra_end[groupLen] = extra_end;
+        overlapID[groupLen] = currentID;
+        y_startGroup[groupLen] = y_start;
+        x_string = r_string + x_start;
+        groupLen++;
         
-        
-        ///这个长度足够，可以用来一起比
-        if (Window_Len == o_len)
+
+        if (groupLen == GROUP_SIZE)
         {
+            Reserve_Banded_BPM_4_SSE_only(dumy->overlap_region_group[0], dumy->overlap_region_group[1], 
+            dumy->overlap_region_group[2], dumy->overlap_region_group[3], Window_Len, x_string, WINDOW,
+                return_sites, return_sites_error, THRESHOLD, dumy->Peq_SSE);
+            groupLen = 0;
 
-            recover_UC_Read_sub_region(dumy->overlap_region_group[groupLen], y_start, o_len, overlap_list->list[currentID].y_pos_strand, 
-                R_INF, overlap_list->list[currentID].y_id);
             
-            overlapID[groupLen] = currentID;
-            y_startGroup[groupLen] = y_start;
-
-            x_string = r_string + x_start;
-
-            groupLen++;
-            if (groupLen == GROUP_SIZE)
+            if (return_sites_error[0]!=(unsigned int)-1)
             {
-            
-                Reserve_Banded_BPM_4_SSE_only(dumy->overlap_region_group[0], dumy->overlap_region_group[1], 
-                dumy->overlap_region_group[2], dumy->overlap_region_group[3], o_len, x_string, x_len,
-			        return_sites, return_sites_error, THRESHOLD, dumy->Peq_SSE);
-                groupLen = 0;
+                overlap_list->list[overlapID[0]].align_length += x_len;
 
-                
-                if (return_sites_error[0]!=(unsigned int)-1)
-                {
-                    overlap_list->list[overlapID[0]].align_length += x_len;
-
-                    append_window_list(&overlap_list->list[overlapID[0]], window_start, window_end, 
-                                y_startGroup[0], y_startGroup[0] + return_sites[0], (int)return_sites_error[0]);
-                }
-                else
-                {
-                    append_window_list(&overlap_list->list[overlapID[0]], window_start, window_end, y_startGroup[0], -1, -1);
-                }
-                
-
-                if (return_sites_error[1]!=(unsigned int)-1)
-                {
-                    overlap_list->list[overlapID[1]].align_length += x_len;
-
-                    append_window_list(&overlap_list->list[overlapID[1]], window_start, window_end, 
-                                y_startGroup[1], y_startGroup[1] + return_sites[1], (int)return_sites_error[1]);
-                }
-                else
-                {
-                    append_window_list(&overlap_list->list[overlapID[1]], window_start, window_end, y_startGroup[1], -1, -1);
-                }
-                
-
-                if (return_sites_error[2]!=(unsigned int)-1)
-                {
-                    overlap_list->list[overlapID[2]].align_length += x_len;
-
-                    append_window_list(&overlap_list->list[overlapID[2]], window_start, window_end, 
-                                y_startGroup[2], y_startGroup[2] + return_sites[2], (int)return_sites_error[2]);
-                }
-                else
-                {
-                    append_window_list(&overlap_list->list[overlapID[2]], window_start, window_end, y_startGroup[2], -1, -1);
-                }
-                
-
-                if (return_sites_error[3]!=(unsigned int)-1)
-                {
-                    overlap_list->list[overlapID[3]].align_length += x_len;
-
-                    append_window_list(&overlap_list->list[overlapID[3]], window_start, window_end, 
-                                y_startGroup[3], y_startGroup[3] + return_sites[3], (int)return_sites_error[3]);
-                }
-                else
-                {
-                    append_window_list(&overlap_list->list[overlapID[3]], window_start, window_end, y_startGroup[3], -1, -1);
-                }
-                
-
-                
-               
-               /************需要注释掉**********/
-            //    for (size_t ijk = 0; ijk < GROUP_SIZE; ijk++)
-            //    {
-            //        if (return_sites_error[ijk]!=(unsigned int)-1)
-            //         {
-            //             total_match++;
-            //         }
-            //         else
-            //         {
-            //             total_unmatch++;
-            //         }
-            //         test_edit_distance_by_edlib(x_string, dumy->overlap_region_group[ijk], x_len, o_len, THRESHOLD, 
-            //         return_sites_error[ijk], &total_mis);
-            //    }
-                /************需要注释掉********* */
-                
-            }
-
-        }
-        else  ///这个不够，只能单个比 ///不够的地方要置N
-        {
-
-            recover_UC_Read_sub_region(dumy->overlap_region, y_start, o_len, overlap_list->list[currentID].y_pos_strand, 
-                R_INF, overlap_list->list[currentID].y_id);
-
-            x_string = r_string + x_start;
-            y_string = dumy->overlap_region;
-
-            ///o_len < Window_Len, 说明y的长度不够，需要在y后面补N
-            memset (y_string + o_len, 0, Window_Len - o_len);
-            end_site = Reserve_Banded_BPM(y_string, o_len, x_string, x_len, THRESHOLD, &error);
-
-            
-
-
-            if (error!=(unsigned int)-1)
-            {
-                overlap_list->list[currentID].align_length += x_len;
-                append_window_list(&overlap_list->list[currentID], window_start, window_end, y_start, y_start + end_site, 
-                (int)error);
+                append_window_list(&overlap_list->list[overlapID[0]], window_start, window_end, 
+                            y_startGroup[0], y_startGroup[0] + return_sites[0], (int)return_sites_error[0],
+                            y_extra_begin[0], y_extra_end[0]);
             }
             else
             {
-                append_window_list(&overlap_list->list[currentID], window_start, window_end, y_start, -1, -1);
+                append_window_list(&overlap_list->list[overlapID[0]], window_start, window_end, y_startGroup[0], -1, -1,
+                y_extra_begin[0], y_extra_end[0]);
             }
             
+
+            if (return_sites_error[1]!=(unsigned int)-1)
+            {
+                overlap_list->list[overlapID[1]].align_length += x_len;
+
+                append_window_list(&overlap_list->list[overlapID[1]], window_start, window_end, 
+                            y_startGroup[1], y_startGroup[1] + return_sites[1], (int)return_sites_error[1],
+                            y_extra_begin[1], y_extra_end[1]);
+            }
+            else
+            {
+                append_window_list(&overlap_list->list[overlapID[1]], window_start, window_end, y_startGroup[1], -1, -1,
+                y_extra_begin[1], y_extra_end[1]);
+            }
             
-            /************需要注释掉**********/
-            // if (error!=(unsigned int)-1)
-            // {
-            //     total_match++;
-            // }
-            // else
-            // {
-            //     total_unmatch++;
-            // }
-            // test_edit_distance_by_edlib(x_string, y_string, x_len, o_len, THRESHOLD, error, &total_mis);
-            /************需要注释掉********* */
+
+            if (return_sites_error[2]!=(unsigned int)-1)
+            {
+                overlap_list->list[overlapID[2]].align_length += x_len;
+
+                append_window_list(&overlap_list->list[overlapID[2]], window_start, window_end, 
+                            y_startGroup[2], y_startGroup[2] + return_sites[2], (int)return_sites_error[2],
+                            y_extra_begin[2], y_extra_end[2]);
+            }
+            else
+            {
+                append_window_list(&overlap_list->list[overlapID[2]], window_start, window_end, y_startGroup[2], -1, -1,
+                y_extra_begin[2], y_extra_end[2]);
+            }
+            
+
+            if (return_sites_error[3]!=(unsigned int)-1)
+            {
+                overlap_list->list[overlapID[3]].align_length += x_len;
+
+                append_window_list(&overlap_list->list[overlapID[3]], window_start, window_end, 
+                            y_startGroup[3], y_startGroup[3] + return_sites[3], (int)return_sites_error[3],
+                            y_extra_begin[3], y_extra_end[3]);
+            }
+            else
+            {
+                append_window_list(&overlap_list->list[overlapID[3]], window_start, window_end, y_startGroup[3], -1, -1,
+                y_extra_begin[3], y_extra_end[3]);
+            }    
         }
     }
 
 
+    
 
     if (groupLen == 1)
     {
-        end_site = Reserve_Banded_BPM(dumy->overlap_region_group[0], o_len, x_string, x_len, THRESHOLD, &error);
+        end_site = Reserve_Banded_BPM(dumy->overlap_region_group[0], Window_Len, x_string, WINDOW, THRESHOLD, &error);
 
         
 
@@ -654,30 +617,19 @@ char* r_string)
             overlap_list->list[overlapID[0]].align_length += x_len;
 
             append_window_list(&overlap_list->list[overlapID[0]], window_start, window_end, 
-                                y_startGroup[0], y_startGroup[0] + end_site, (int)error);
+                                y_startGroup[0], y_startGroup[0] + end_site, (int)error,
+                                y_extra_begin[0], y_extra_end[0]);
         }
         else
         {
-            append_window_list(&overlap_list->list[overlapID[0]], window_start, window_end, y_startGroup[0], -1, -1);
+            append_window_list(&overlap_list->list[overlapID[0]], window_start, window_end, y_startGroup[0], -1, -1,
+            y_extra_begin[0], y_extra_end[0]);
         }
-        
-
-        /************需要注释掉**********/
-        // if (error!=(unsigned int)-1)
-        // {
-        //     total_match++;
-        // }
-        // else
-        // {
-        //     total_unmatch++;
-        // }
-        // test_edit_distance_by_edlib(x_string, dumy->overlap_region_group[0], x_len, o_len, THRESHOLD, error, &total_mis);
-        /************需要注释掉********* */
     }
     else if (groupLen > 1)
     {
         Reserve_Banded_BPM_4_SSE_only(dumy->overlap_region_group[0], dumy->overlap_region_group[1], 
-                dumy->overlap_region_group[2], dumy->overlap_region_group[3], o_len, x_string, x_len,
+                dumy->overlap_region_group[2], dumy->overlap_region_group[3], Window_Len, x_string, WINDOW,
 			        return_sites, return_sites_error, THRESHOLD, dumy->Peq_SSE);
 
         for (i = 0; i < groupLen; i++)
@@ -686,35 +638,20 @@ char* r_string)
             {
                 overlap_list->list[overlapID[i]].align_length += x_len;
                 append_window_list(&overlap_list->list[overlapID[i]], window_start, window_end, 
-                                y_startGroup[i], y_startGroup[i] + return_sites[i], (int)return_sites_error[i]);
+                                y_startGroup[i], y_startGroup[i] + return_sites[i], (int)return_sites_error[i],
+                                y_extra_begin[i], y_extra_end[i]);
             }
             else
             {
-                append_window_list(&overlap_list->list[overlapID[i]], window_start, window_end, y_startGroup[i], -1, -1);
+                append_window_list(&overlap_list->list[overlapID[i]], window_start, window_end, y_startGroup[i], -1, -1,
+                y_extra_begin[i], y_extra_end[i]);
             }
             
         }
 
-
-        /************需要注释掉**********/
-        // for (size_t ijk = 0; ijk < groupLen; ijk++)
-        // {
-        //     if (return_sites_error[ijk]!=(unsigned int)-1)
-        //     {
-        //         total_match++;
-        //     }
-        //     else
-        //     {
-        //         total_unmatch++;
-        //     }
-        //     test_edit_distance_by_edlib(x_string, dumy->overlap_region_group[ijk], x_len, o_len, THRESHOLD, 
-        //     return_sites_error[ijk], &total_mis);
-        // }
-        /************需要注释掉**********/
         groupLen = 0;
     }
     
-    
 
     long long reverse_i = dumy->size - 1;
     int threshold;
@@ -722,6 +659,7 @@ char* r_string)
     ///这些是整个window被部分覆盖的
     for (i = 0; i < dumy->lengthNT; i++)
     {
+        extra_begin = extra_end = 0;
         currentID = dumy->overlapID[reverse_i--];
         x_start = MAX(window_start, overlap_list->list[currentID].x_pos_s);
         x_end = MIN(window_end, overlap_list->list[currentID].x_pos_e);
@@ -730,104 +668,38 @@ char* r_string)
         x_len = x_end - x_start + 1;
         threshold = x_len * THRESHOLD_RATE;
 
-        // /************需要注释掉**********/
-        // if (x_len <= 0)
-        // {
-        //     fprintf(stderr, "ERROR\n");
-        // }
-        // /************需要注释掉**********/
-
         ///y上的相对位置
         y_start = (x_start - overlap_list->list[currentID].x_pos_s) + overlap_list->list[currentID].y_pos_s;
 
-        // /************需要注释掉**********/
-        // if(x_start < overlap_list->list[currentID].x_pos_s)
-        // {
-        //     fprintf(stderr, "ERROR\n");
-        // }
-        // /************需要注释掉**********/
-
-
-        ///y上的起始
-        ///y_start = y_start - THRESHOLD;
-        y_start = y_start - threshold;
-        if (y_start < 0)
-        {
-            y_start = 0;
-        }
         
-        ///当前y的长度
-        currentIDLen = Get_READ_LENGTH((*R_INF), overlap_list->list[currentID].y_id);
-
-        
-        ///不能超过y的剩余长度
-        ///Window_Len = x_len + (THRESHOLD << 1);
         Window_Len = x_len + (threshold << 1);
-        o_len = MIN(Window_Len, currentIDLen - y_start);
-
-
-        // /************需要注释掉**********/
-        // if(o_len < x_len)
-        // {
-        //     fprintf(stderr, "ERROR\n");
-        // }
-        // /************需要注释掉**********/
+        determine_overlap_region(threshold, y_start, overlap_list->list[currentID].y_id, Window_Len, R_INF,
+        &extra_begin, &extra_end, &y_start, &o_len);
 
+        fill_subregion(dumy->overlap_region, y_start, o_len, overlap_list->list[currentID].y_pos_strand, 
+        R_INF, overlap_list->list[currentID].y_id, extra_begin, extra_end);
         
-        recover_UC_Read_sub_region(dumy->overlap_region, y_start, o_len, overlap_list->list[currentID].y_pos_strand, 
-                R_INF, overlap_list->list[currentID].y_id);
 
         x_string = r_string + x_start;
         y_string = dumy->overlap_region;
 
-        ///不够的地方要置N
-        if (Window_Len != o_len)
-        {
-            ///o_len < Window_Len, 说明y的长度不够，需要在y后面补N
-            memset (y_string + o_len, 0, Window_Len - o_len);
-        }
-        ///end_site = Reserve_Banded_BPM(y_string, o_len, x_string, x_len, THRESHOLD, &error);
-        end_site = Reserve_Banded_BPM(y_string, o_len, x_string, x_len, threshold, &error);
+
+        end_site = Reserve_Banded_BPM(y_string, Window_Len, x_string, x_len, threshold, &error);
 
 
         if (error!=(unsigned int)-1)
         {
             overlap_list->list[currentID].align_length += x_len;
-            append_window_list(&overlap_list->list[currentID], x_start, x_end, y_start, y_start + end_site, (int)error);
+            append_window_list(&overlap_list->list[currentID], x_start, x_end, y_start, y_start + end_site, (int)error,
+            extra_begin, extra_end);
         }
         else
         {
-            append_window_list(&overlap_list->list[currentID], x_start, x_end, y_start, -1, -1);
+            append_window_list(&overlap_list->list[currentID], x_start, x_end, y_start, -1, -1,
+            extra_begin, extra_end);
         }
-
-
-
-
-        /************需要注释掉********* */
-        // if (error!=(unsigned int)-1)
-        // {
-        //     total_match++;
-        // }
-        // else
-        // {
-        //     total_unmatch++;
-        // }
-        // //test_edit_distance_by_edlib(x_string, y_string, x_len, o_len, THRESHOLD, error, &total_mis);
-        // test_edit_distance_by_edlib(x_string, y_string, x_len, o_len, threshold, error, &total_mis);
-        /************需要注释掉********* */
     }
-    
 
-
-
-    /************需要注释掉********* */
-    // pthread_mutex_lock(&debug_statistics);
-    // T_total_match = T_total_match + total_match;
-    // T_total_unmatch = T_total_unmatch + total_unmatch;
-    // T_total_mis = T_total_mis + total_mis;
-	// pthread_mutex_unlock(&debug_statistics);
-    /************需要注释掉********* */
-        
 }
 
 void debug_stats(overlap_region_alloc* overlap_list, All_reads* R_INF, 
@@ -933,7 +805,7 @@ void debug_stats(overlap_region_alloc* overlap_list, All_reads* R_INF,
 
 
 inline void generate_cigar(
-	char* path, int path_length, window_list* result)
+	char* path, int path_length, window_list* result, int* start, int* end)
 {
 
     if (result->error == 0)
@@ -952,6 +824,34 @@ inline void generate_cigar(
     char pre_ciga = 5;
     int pre_ciga_length = 0;
 
+
+    for (i = 0; i < path_length; i++)
+    {
+        if(path[i] == 1)
+        {
+            path[i] = 3;
+            (*end)--;
+        }
+        else
+        {
+            break;
+        }
+    }
+
+    for (i = path_length - 1; i >= 0; i--)
+    {
+        if(path[i] == 1)
+        {
+            path[i] = 3;
+            (*start)++;
+        }
+        else
+        {
+            break;
+        }
+    }
+    
+
     for (i = path_length - 1; i >= 0; i--)
     {
 
@@ -1110,6 +1010,8 @@ inline void recalcate_window(overlap_region_alloc* overlap_list, All_reads* R_IN
     unsigned int error;
     int real_y_start;
     long long overlap_length;
+    int extra_begin, extra_end;
+    long long o_len;
 
 
     ///j负责遍历整个overlap list
@@ -1130,12 +1032,20 @@ inline void recalcate_window(overlap_region_alloc* overlap_list, All_reads* R_IN
             ///找到第一个匹配的window
             if(overlap_list->list[j].w_list[i].y_end != -1)
             {
-                total_y_start = overlap_list->list[j].w_list[i].y_end + 1;
+                ///note!!! need notification
+                ///total_y_start = overlap_list->list[j].w_list[i].y_end + 1;
+                ///this is the actual end postion in ystring
+                total_y_start = overlap_list->list[j].w_list[i].y_end - overlap_list->list[j].w_list[i].extra_begin + 1;
+
+
+
                 ///k遍历匹配window右侧所有不匹配的window
                 ///如果i匹配，则k从i+1开始
                 ///知道第一个匹配的window结束
                 for (k = i + 1; k < overlap_list->list[j].w_list_length && overlap_list->list[j].w_list[k].y_end == -1; k++)
-                {   
+                {
+                    extra_begin = extra_end = 0;
+
                     ///y_start有可能大于y_readLen
                     ///这多发于最后一个window长度仅为几，而前面一个window的结束位置也超过了y_readLen-1
                     ///这个时候做动态规划会给超过的部分补N
@@ -1144,40 +1054,31 @@ inline void recalcate_window(overlap_region_alloc* overlap_list, All_reads* R_IN
                         break;
                     }
                     
+                    ///there is no problem for x
                     x_start = overlap_list->list[j].w_list[k].x_start;
                     x_end = overlap_list->list[j].w_list[k].x_end;
                     x_len = x_end - x_start + 1;
-                    ///if(x_len == )
                     threshold = x_len * THRESHOLD_RATE;
-                    y_start = total_y_start - threshold;
-                    if (y_start < 0)
-                    {
-                        y_start = 0;
-                    }
+
+
+                    y_start = total_y_start;
                     Window_Len = x_len + (threshold << 1);
-                    ///y_start有可能大于y_readLen
-                    ///这多发于最后一个window长度仅为几，而前面一个window的结束位置也超过了y_readLen-1
-                    ///这个时候做动态规划会给超过的部分补N
-                    y_len = MIN(Window_Len, y_readLen - y_start);
-                    ///这说明已经到y的结尾了
-                    if (y_len < x_len)
+                    determine_overlap_region(threshold, y_start, y_id, Window_Len, R_INF, 
+                    &extra_begin, &extra_end, &y_start, &o_len);
+
+                    if(o_len + threshold < x_len)
                     {
                         break;
                     }
                     
-                    
+                    fill_subregion(dumy->overlap_region, y_start, o_len, y_strand, 
+                    R_INF, y_id, extra_begin, extra_end);
 
-                    recover_UC_Read_sub_region(dumy->overlap_region, y_start, y_len, y_strand, R_INF, y_id);
-                    
                     x_string = g_read->seq + x_start;
                     y_string = dumy->overlap_region;
-                    if (Window_Len != y_len)
-                    {
-                        ///o_len < Window_Len, 说明y的长度不够，需要在y后面补N
-                        memset (y_string + y_len, 0, Window_Len - y_len);
-                    }
-                    
-                    end_site = Reserve_Banded_BPM(y_string, y_len, x_string, x_len, threshold, &error);
+
+                    ///note!!! need notification
+                    end_site = Reserve_Banded_BPM(y_string, Window_Len, x_string, x_len, threshold, &error);
                     
                     ///error等于-1说明没匹配
                     if (error!=(unsigned int)-1)
@@ -1186,14 +1087,18 @@ inline void recalcate_window(overlap_region_alloc* overlap_list, All_reads* R_IN
                         overlap_list->list[j].w_list[k].y_start = y_start;
                         overlap_list->list[j].w_list[k].y_end = y_start + end_site;
                         overlap_list->list[j].w_list[k].error = (int)error;
+                        ///note!!! need notification
+                        overlap_list->list[j].w_list[k].extra_begin = extra_begin;
+                        overlap_list->list[j].w_list[k].extra_end = extra_end;
                         overlap_list->list[j].align_length += x_len;
                     }
                     else
                     {
                         break;
                     }
-                    
-                    total_y_start = y_start + end_site + 1;
+                    ///note!!! need notification
+                    ///total_y_start = y_start + end_site + 1;
+                    total_y_start = y_start + end_site - extra_begin + 1;
 
                 }
                 
@@ -1222,25 +1127,26 @@ inline void recalcate_window(overlap_region_alloc* overlap_list, All_reads* R_IN
                 ///如果没有，就需要重新计算
                 if(overlap_list->list[j].w_list[i].cigar.length == -1)
                 {
+                    ///there is no problem for x
                     x_start = overlap_list->list[j].w_list[i].x_start;
                     x_end = overlap_list->list[j].w_list[i].x_end;
                     x_len = x_end - x_start + 1;
                     threshold = x_len * THRESHOLD_RATE;
-                    y_start = overlap_list->list[j].w_list[i].y_start;
                     Window_Len = x_len + (threshold << 1);
-                    y_len = MIN(Window_Len, y_readLen - y_start);
 
-                    recover_UC_Read_sub_region(dumy->overlap_region, y_start, y_len, y_strand, R_INF, y_id);
 
+                    ///y_start is the real y_start
+                    y_start = overlap_list->list[j].w_list[i].y_start;
+                    extra_begin = overlap_list->list[j].w_list[i].extra_begin;
+                    extra_end = overlap_list->list[j].w_list[i].extra_end;
+                    o_len = Window_Len - extra_end - extra_begin;
+                    fill_subregion(dumy->overlap_region, y_start, o_len, y_strand, 
+                    R_INF, y_id, extra_begin, extra_end);
                     x_string = g_read->seq + x_start;
                     y_string = dumy->overlap_region;
-                    if (Window_Len != y_len)
-                    {
-                        ///o_len < Window_Len, 说明y的长度不够，需要在y后面补N
-                        memset (y_string + y_len, 0, Window_Len - y_len);
-                    }
 
-                    end_site = Reserve_Banded_BPM_PATH(y_string, y_len, x_string, x_len, threshold, &error, &real_y_start,
+                    ///note!!! need notification
+                    end_site = Reserve_Banded_BPM_PATH(y_string, Window_Len, x_string, x_len, threshold, &error, &real_y_start,
                     &(dumy->path_length), dumy->matrix_bit, dumy->path, 
                     overlap_list->list[j].w_list[i].error, overlap_list->list[j].w_list[i].y_end - y_start);
 
@@ -1250,10 +1156,34 @@ inline void recalcate_window(overlap_region_alloc* overlap_list, All_reads* R_IN
                     ///到这里y_start已经被正确计算出来了
                     if (error != (unsigned int)-1)
                     {
-                        real_y_start = y_start + real_y_start;
-                        overlap_list->list[j].w_list[i].y_start = real_y_start;
-                        generate_cigar(dumy->path, dumy->path_length, &(overlap_list->list[j].w_list[i]));                                
+                        
+
+
+                         
+                        generate_cigar(dumy->path, dumy->path_length, &(overlap_list->list[j].w_list[i]),
+                        &real_y_start, &end_site);  
+
+                        if(real_y_start < extra_begin || end_site >= Window_Len - extra_end)
+                        {
+                            fprintf(stderr, "\nreal_y_start: %d, extra_begin: %d\n", 
+                            real_y_start, extra_begin);
+
+                            fprintf(stderr, "end_site: %d, Window_Len: %d, extra_end: %d\n", 
+                            end_site, Window_Len, extra_end);                            
+                        }
+                        
+                        
+                        ///note!!! need notification
+                        ///real_y_start = y_start + real_y_start;
+                        real_y_start = y_start + real_y_start - extra_begin;
+                        overlap_list->list[j].w_list[i].y_start = real_y_start;  
+                        overlap_list->list[j].w_list[i].y_end = y_start + end_site;                         
                     }
+                    else
+                    {
+                        fprintf(stderr, "error\n");
+                    }
+                    
                         
                 }
                 else
@@ -1271,59 +1201,69 @@ inline void recalcate_window(overlap_region_alloc* overlap_list, All_reads* R_IN
                 ///因为i!=0，所以k的大小不用担心
                 for (k = i - 1; k >= 0 && overlap_list->list[j].w_list[k].y_end == -1; k--)
                 {  
+                    ///there is no problem in x
                     x_start = overlap_list->list[j].w_list[k].x_start;
                     x_end = overlap_list->list[j].w_list[k].x_end;
                     x_len = x_end - x_start + 1;
                     threshold = x_len * THRESHOLD_RATE;
-                    ///这个和上面不同，先求y_end
-                    y_end = total_y_end + threshold;
-                    ///y_end不能大于y的总长度
-                    if(y_end >= y_readLen)
-                    {
-                        y_end = y_readLen - 1;
-                    }
                     Window_Len = x_len + (threshold << 1);
-                    y_len = MIN(Window_Len, y_end + 1);
-                    ///这说明已经到y的开始了，因为是倒着算的，没必要接着算了
-                    if (y_len < x_len)
+
+                    if(total_y_end <= 0)
                     {
                         break;
                     }
-                    y_start = y_end - y_len + 1;
-                    recover_UC_Read_sub_region(dumy->overlap_region, y_start, y_len, y_strand, R_INF, y_id);
+
+                    ///y_start may less than 0
+                    y_start = total_y_end - x_len + 1;
+                    determine_overlap_region(threshold, y_start, y_id, Window_Len, R_INF,
+                    &extra_begin, &extra_end, &y_start, &o_len);
+
+                    if(o_len + threshold < x_len)
+                    {
+                        break;
+                    }
+
+                    fill_subregion(dumy->overlap_region, y_start, o_len, y_strand, 
+                    R_INF, y_id, extra_begin, extra_end);
                     x_string = g_read->seq + x_start;
                     y_string = dumy->overlap_region;
-                    
-                    ///不要在前面补补，补了有可能触发剪纸
-                    ///还是像上面一样在后面随便补点就好了
-                    if (Window_Len != y_len)
-                    {
-                        ///o_len < Window_Len, 说明y的长度不够，需要在y后面补N
-                        ///memset (y_string, 0, Window_Len - y_len);
-                        memset (y_string + y_len, 0, Window_Len - y_len);
-                    }
-                    
 
-                    end_site = Reserve_Banded_BPM_PATH(y_string, y_len, x_string, x_len, threshold, &error, &real_y_start,
+                    ///note!!! need notification
+                    end_site = Reserve_Banded_BPM_PATH(y_string, Window_Len, x_string, x_len, threshold, &error, &real_y_start,
                     &(dumy->path_length), dumy->matrix_bit, dumy->path, -1, -1);
 
+
                     ///error等于-1说明没匹配
                     if (error!=(unsigned int)-1)
-                    {
-                        
-                        ///overlap_list->list[j].w_list[k].y_pre_start = overlap_list->list[j].w_list[k].y_start;
-                        overlap_list->list[j].w_list[k].y_start = y_start + real_y_start;
+                    { 
+
+                        generate_cigar(dumy->path, dumy->path_length, &(overlap_list->list[j].w_list[k]),
+                        &real_y_start, &end_site);   
+
+                        if(real_y_start < extra_begin || end_site >= Window_Len - extra_end)
+                        {
+                            fprintf(stderr, "\nreal_y_start: %d, extra_begin: %d\n", 
+                            real_y_start, extra_begin);
+
+                            fprintf(stderr, "end_site: %d, Window_Len: %d, extra_end: %d\n", 
+                            end_site, Window_Len, extra_end);                            
+                        }
+
+                        ///y_start has no shift, but y_end has shift           
+                        overlap_list->list[j].w_list[k].y_start = y_start + real_y_start - extra_begin;
                         overlap_list->list[j].w_list[k].y_end = y_start + end_site;
                         overlap_list->list[j].w_list[k].error = error;
                         overlap_list->list[j].align_length += x_len;
-                        generate_cigar(dumy->path, dumy->path_length, &(overlap_list->list[j].w_list[k]));
+                        ///note!!! need notification
+                        overlap_list->list[j].w_list[k].extra_begin = extra_begin;
+                        overlap_list->list[j].w_list[k].extra_end = extra_end;
                     }
                     else
                     {
                         break;
                     }
 
-                    total_y_end = y_start + real_y_start - 1;
+                    total_y_end = y_start + real_y_start - 1 - extra_begin;
                 }
             
             }
@@ -1344,6 +1284,7 @@ inline void recalcate_window(overlap_region_alloc* overlap_list, All_reads* R_IN
         y_readLen = Get_READ_LENGTH((*R_INF), y_id);
         overlap_length = overlap_list->list[j].x_pos_e - overlap_list->list[j].x_pos_s + 1;
 
+        ///only calculate cigar for high quality overlaps
         if (overlap_length * OVERLAP_THRESHOLD <=  overlap_list->list[j].align_length)
         {
             for (i = 0; i < overlap_list->list[j].w_list_length; i++)
@@ -1351,38 +1292,72 @@ inline void recalcate_window(overlap_region_alloc* overlap_list, All_reads* R_IN
                 ///判断cigar是否被计算
                 ///没被计算过就重算
                 ///第一个条件是判断这个窗口是否匹配
-                if(overlap_list->list[j].w_list[i].y_end != -1 && overlap_list->list[j].w_list[i].cigar.length == -1)
+                if(overlap_list->list[j].w_list[i].y_end != -1)
                 {
-                    x_start = overlap_list->list[j].w_list[i].x_start;
-                    x_end = overlap_list->list[j].w_list[i].x_end;
-                    x_len = x_end - x_start + 1;
-                    threshold = x_len * THRESHOLD_RATE;
-                    y_start = overlap_list->list[j].w_list[i].y_start;
-                    Window_Len = x_len + (threshold << 1);
-                    y_len = MIN(Window_Len, y_readLen - y_start);
-
-                    recover_UC_Read_sub_region(dumy->overlap_region, y_start, y_len, y_strand, R_INF, y_id);
-
-                    x_string = g_read->seq + x_start;
-                    y_string = dumy->overlap_region;
-                    if (Window_Len != y_len)
+                    if(overlap_list->list[j].w_list[i].cigar.length == -1)
                     {
-                        ///o_len < Window_Len, 说明y的长度不够，需要在y后面补N
-                        memset (y_string + y_len, 0, Window_Len - y_len);
+                        ///there is no problem for x
+                        x_start = overlap_list->list[j].w_list[i].x_start;
+                        x_end = overlap_list->list[j].w_list[i].x_end;
+                        x_len = x_end - x_start + 1;
+                        threshold = x_len * THRESHOLD_RATE;
+                        Window_Len = x_len + (threshold << 1);
+
+
+                        ///y_start is the real y_start
+                        y_start = overlap_list->list[j].w_list[i].y_start;
+                        extra_begin = overlap_list->list[j].w_list[i].extra_begin;
+                        extra_end = overlap_list->list[j].w_list[i].extra_end;
+                        o_len = Window_Len - extra_end - extra_begin;
+                        fill_subregion(dumy->overlap_region, y_start, o_len, y_strand, 
+                        R_INF, y_id, extra_begin, extra_end);
+                        x_string = g_read->seq + x_start;
+                        y_string = dumy->overlap_region;
+
+
+                        ///note!!! need notification
+                        end_site = Reserve_Banded_BPM_PATH(y_string, Window_Len, x_string, x_len, threshold, &error, &real_y_start,
+                        &(dumy->path_length), dumy->matrix_bit, dumy->path, 
+                        overlap_list->list[j].w_list[i].error, overlap_list->list[j].w_list[i].y_end - y_start);
+
+
+                        ///到这里y_start已经被正确计算出来了
+                        if (error != (unsigned int)-1)
+                        {
+
+                            generate_cigar(dumy->path, dumy->path_length, &(overlap_list->list[j].w_list[i]),
+                            &real_y_start, &end_site);    
+
+
+                            if(real_y_start < extra_begin || end_site >= Window_Len - extra_end)
+                            {
+                                fprintf(stderr, "\nreal_y_start: %d, extra_begin: %d\n", 
+                                real_y_start, extra_begin);
+
+                                fprintf(stderr, "end_site: %d, Window_Len: %d, extra_end: %d\n", 
+                                end_site, Window_Len, extra_end);                            
+                            }
+
+                            
+
+                            ///note!!! need notification
+                            ///real_y_start = y_start + real_y_start;
+                            real_y_start = y_start + real_y_start - extra_begin;
+                            overlap_list->list[j].w_list[i].y_start = real_y_start;  
+                            overlap_list->list[j].w_list[i].y_end = y_start + end_site - extra_begin;                              
+                        }
+                        else
+                        {
+                            fprintf(stderr, "error\n");
+                        }
                     }
-
-                    end_site = Reserve_Banded_BPM_PATH(y_string, y_len, x_string, x_len, threshold, &error, &real_y_start,
-                    &(dumy->path_length), dumy->matrix_bit, dumy->path, 
-                    overlap_list->list[j].w_list[i].error, overlap_list->list[j].w_list[i].y_end - y_start);
-
-                    ///到这里y_start已经被正确计算出来了
-                    if (error != (unsigned int)-1)
+                    else
                     {
-                        real_y_start = y_start + real_y_start;
-                        overlap_list->list[j].w_list[i].y_start = real_y_start;
-                        generate_cigar(dumy->path, dumy->path_length, &(overlap_list->list[j].w_list[i]));                                
+                        overlap_list->list[j].w_list[i].y_end -= overlap_list->list[j].w_list[i].extra_begin;
                     }
                 }
+                
+                
             }
         }
     }
@@ -1391,205 +1366,48 @@ inline void recalcate_window(overlap_region_alloc* overlap_list, All_reads* R_IN
 
 
 
+    /**
+    ///j负责遍历整个overlap list
+    for (j = 0; j < overlap_list->length; j++)
+    {
+        y_id = overlap_list->list[j].y_id;
+        y_strand = overlap_list->list[j].y_pos_strand;
+        y_readLen = Get_READ_LENGTH((*R_INF), y_id);
+        overlap_length = overlap_list->list[j].x_pos_e - overlap_list->list[j].x_pos_s + 1;
+
+        ///only calculate cigar for high quality overlaps
+        if (overlap_length * OVERLAP_THRESHOLD <=  overlap_list->list[j].align_length)
+        {
+            for (i = 0; i < overlap_list->list[j].w_list_length; i++)
+            {
+                if(overlap_list->list[j].w_list[i].y_end != -1)
+                {
+                    ///there is no problem for x
+                    x_start = overlap_list->list[j].w_list[i].x_start;
+                    x_end = overlap_list->list[j].w_list[i].x_end;
+                    x_len = x_end - x_start + 1;
+
+                    x_string = g_read->seq + x_start;
+
+                    y_start = overlap_list->list[j].w_list[i].y_start;
+                    y_end = overlap_list->list[j].w_list[i].y_end;
+                    y_len = y_end - y_start + 1;
+
+                    recover_UC_Read_sub_region(dumy->overlap_region, y_start, y_len, y_strand, R_INF, y_id);
+                    y_string = dumy->overlap_region;
 
 
+                    if(verify_cigar(x_string, x_len, y_string, y_len, &overlap_list->list[j].w_list[i].cigar, 
+                    overlap_list->list[j].w_list[i].error))
+                    {
+                        fprintf(stderr, "j: %d, i: %d, y_id: %d, y_start: %d, y_end: %d\n", j, i, y_id, y_start, y_end);
+                    }
+                }
 
-
-
-
-
-
-
-
-
-
-    
-
-    
-    // for (j = 0; j < overlap_list->length; j++)
-    // {
-        
-    //     y_id = overlap_list->list[j].y_id;
-    //     y_strand = overlap_list->list[j].y_pos_strand;
-    //     y_readLen = Get_READ_LENGTH((*R_INF), y_id);
-
-    //     matches = 0;
-    //     for (i = 0; i < overlap_list->list[j].w_list_length; i++)
-    //     {
-    //         x_start = overlap_list->list[j].w_list[i].x_start;
-    //         x_end = overlap_list->list[j].w_list[i].x_end;
-    //         x_len = x_end - x_start + 1;
-    //         ///if(x_len == )
-    //         threshold = x_len * THRESHOLD_RATE;   
-    //         y_start = overlap_list->list[j].w_list[i].y_start;
-
-
-            
-    //         if (overlap_list->list[j].w_list[i].y_end != -1)
-    //         {
-    //             matches += x_len;
-    //         }
-
-    //         if (i + 1 < overlap_list->list[j].w_list_length)
-    //         {
-    //             if (
-    //                 overlap_list->list[j].w_list[i + 1].x_end <= overlap_list->list[j].w_list[i].x_end
-    //                 ||
-    //                 overlap_list->list[j].w_list[i + 1].x_start <= overlap_list->list[j].w_list[i].x_start
-    //                 )
-    //             {
-    //                 fprintf(stderr, "ERROR 1\n");
-    //             }
-                
-    //         }
-
-
-    //         ///当起始位置还没求出来的时候
-    //         if(overlap_list->list[j].w_list[i].cigar.length == -1)
-    //         {
-                
-    //             Window_Len = x_len + (threshold << 1);
-    //             y_len = MIN(Window_Len, y_readLen - y_start);
-
-    //             recover_UC_Read_sub_region(dumy->overlap_region, y_start, y_len, y_strand, R_INF, y_id);
-                    
-    //             x_string = g_read->seq + x_start;
-    //             y_string = dumy->overlap_region;
-    //             if (Window_Len != y_len)
-    //             {
-    //                 ///o_len < Window_Len, 说明y的长度不够，需要在y后面补N
-    //                 memset (y_string + y_len, 0, Window_Len - y_len);
-    //             }
-
-    //             ///匹配
-    //             if (overlap_list->list[j].w_list[i].y_end != -1)
-    //             {
-    //                 ///fprintf(stderr, "hahah\n");
-    //                 end_site = Reserve_Banded_BPM(y_string, y_len, x_string, x_len, threshold, &error);
-
-    //                 if (error == (unsigned int)-1)
-    //                 {
-    //                     fprintf(stderr, "1 ERROR, y_start: %u, y_end: %u\n", 
-    //                     overlap_list->list[j].w_list[i].y_start, overlap_list->list[j].w_list[i].y_end);
-    //                 }
-                    
-    //                 if (end_site + y_start != overlap_list->list[j].w_list[i].y_end)
-    //                 {
-    //                     fprintf(stderr, "ERROR, y_end: %u, end_site: %u, y_start: %u, +: %u\n", 
-    //                     overlap_list->list[j].w_list[i].y_end, end_site, y_start, end_site + y_start);
-    //                 }
-
-    //                 int old_error = overlap_list->list[j].w_list[i].error;
-
-    //                 if (error != old_error)
-    //                 {
-                        
-    //                     fprintf(stderr, "ERROR new error: %d, old error: %d\n", error, overlap_list->list[j].w_list[i].error);
-    //                 }
-                    
-                    
-
-    //             }
-    //             else
-    //             {
-    //                 end_site = Reserve_Banded_BPM(y_string, y_len, x_string, x_len, threshold, &error);
-
-    //                 if (error != (unsigned int)-1)
-    //                 {
-    //                     fprintf(stderr, "2 ERROR\n");
-    //                 }
-    //             }
-                
-
-    //         }
-    //         else
-    //         {
-                
-
-    //             y_start = overlap_list->list[j].w_list[i].y_start;
-    //             ///此时y_start是真正匹配的起始位置，所以要减去threshold
-    //             y_start = y_start - threshold;
-    //             if(y_start < 0)
-    //             {
-    //                 y_start = 0;
-    //             }
-
-    //             Window_Len = x_len + (threshold << 1);
-    //             y_len = MIN(Window_Len, y_readLen - y_start);
-
-    //             recover_UC_Read_sub_region(dumy->overlap_region, y_start, y_len, y_strand, R_INF, y_id);
-                    
-    //             x_string = g_read->seq + x_start;
-    //             y_string = dumy->overlap_region;
-
-    //             if (Window_Len != y_len)
-    //             {
-    //                 ///o_len < Window_Len, 说明y的长度不够，需要在y后面补N
-    //                 memset (y_string + y_len, 0, Window_Len - y_len);
-    //             }
-
-
-                
-    //             end_site = Reserve_Banded_BPM_PATH(y_string, y_len, x_string, x_len, threshold, &error, &real_y_start,
-    //                     &(dumy->path_length), dumy->matrix_bit, dumy->path, 
-    //                     overlap_list->list[j].w_list[i].error, overlap_list->list[j].w_list[i].y_end - y_start);
-                        
-    //             ///end_site = Reserve_Banded_BPM(y_string, y_len, x_string, x_len, threshold, &error);
-                
-
-    //             if (error == (unsigned int)-1)
-    //             {
-    //                 fprintf(stderr, "3 ERROR, y_start: %d, y_end: %d, x_len: %d, pre_error: %d, y_strand: %u\n", 
-    //                 overlap_list->list[j].w_list[i].y_start, overlap_list->list[j].w_list[i].y_end, x_len, overlap_list->list[j].w_list[i].error, y_strand);
-
-    //                 fprintf(stderr, "Window_Len: %d, y_len: %d, y_readLen: %d\n", 
-    //                 Window_Len, y_len, y_readLen);
-
-    //                 print_string(g_read->seq+overlap_list->list[j].w_list[i].x_start, 
-    //                 overlap_list->list[j].w_list[i].x_end-overlap_list->list[j].w_list[i].x_start+1);
-
-    //                 recover_UC_Read_sub_region(dumy->overlap_region, overlap_list->list[j].w_list[i].y_start, 
-    //                 overlap_list->list[j].w_list[i].y_end-overlap_list->list[j].w_list[i].y_start + 1, y_strand, R_INF, y_id);
-
-    //                 print_string(dumy->overlap_region, 
-    //                 overlap_list->list[j].w_list[i].y_end-overlap_list->list[j].w_list[i].y_start + 1);
-    //             }
-
-    //             ////检验cigar是否正确
-    //             x_start = overlap_list->list[j].w_list[i].x_start;
-    //             x_end = overlap_list->list[j].w_list[i].x_end;
-    //             x_len = x_end - x_start + 1;
-    //             x_string = g_read->seq + x_start;
-
-    //             y_start = overlap_list->list[j].w_list[i].y_start;
-    //             y_end = overlap_list->list[j].w_list[i].y_end;
-    //             y_len = y_end - y_start + 1;
-    //             recover_UC_Read_sub_region(dumy->overlap_region, y_start, y_len, y_strand, R_INF, y_id);
-    //             y_string = dumy->overlap_region;
-                
-                
-    //             if(verify_cigar(x_string, x_len, y_string, y_len, &overlap_list->list[j].w_list[i].cigar, 
-    //             overlap_list->list[j].w_list[i].error))
-    //             {
-    //                 fprintf(stderr, "j: %d, i: %d, y_id: %d, y_start: %d, y_end: %d\n", j, i, y_id, y_start, y_end);
-    //             }
-                
-                
-                
-            
-                
-
-    //         }
-            
-    //     }
-
-    //     if (matches != overlap_list->list[j].align_length)
-    //     {
-    //         fprintf(stderr, "ERROR 2: matches: %d, align_length: %d\n", 
-    //                                 matches, overlap_list->list[j].align_length);
-    //     }
-    // }
-    
+            }
+        }
+    }
+    **/
     
 }
 
@@ -1654,11 +1472,14 @@ inline void add_segment_to_correct_read(Correct_dumy* dumy, char* segment, long
 
 ///返回下一个backbone节点上的ID
 long long inline add_path_to_correct_read(Graph* backbone, Correct_dumy* dumy, long long currentNodeID, 
-long long type, long long edgeID)
+long long type, long long edgeID, Cigar_record* current_cigar, char* self_string)
 {
     //long long i;
     long long nodeID;
 
+    ///Note: currentNodeID must be a backbone node
+    ///currentNodeID = 0 means a fake node
+    ///currentNodeID = i means self_string[i - 1]
     ///包括匹配和误配两种情况
     if (type == MISMATCH)
     {
@@ -1669,6 +1490,10 @@ long long type, long long edgeID)
             add_base_to_correct_read_directly(dumy, backbone->g_nodes.list[nodeID].base);
             ///match所以dumy->corrected_base不要+1
 
+            ///nodeID = i means self_string[i - 1]
+            ///add_cigar_record(self_string+nodeID-1, 1, current_cigar, 0);
+            add_cigar_record(&(backbone->g_nodes.list[nodeID].base), 1, current_cigar, 0);
+
             /***********需要注释掉********* */
             if (nodeID != currentNodeID + 1)
             {
@@ -1680,12 +1505,26 @@ long long type, long long edgeID)
         }
         else ///这是mismatch的情况
         {
+            
+
             nodeID = backbone->g_nodes.list[currentNodeID].mismatch_edges.list[edgeID].out_node;
             add_base_to_correct_read_directly(dumy, backbone->g_nodes.list[nodeID].base);
             dumy->corrected_base++;
 
+            
+            char merge_base = 0;
+            merge_base = seq_nt6_table[(uint8_t)backbone->g_nodes.list[nodeID].base];
+            merge_base = merge_base << 3;
             ///这种中间节点只有一个元素，所以直接list[0]
             nodeID = backbone->g_nodes.list[nodeID].mismatch_edges.list[0].out_node;
+            merge_base = merge_base | seq_nt6_table[(uint8_t)backbone->g_nodes.list[nodeID].base];
+            add_cigar_record(&merge_base, 1, current_cigar, 1);
+            
+            /**
+            add_cigar_record(&(backbone->g_nodes.list[nodeID].base), 1, current_cigar, 1);
+            nodeID = backbone->g_nodes.list[nodeID].mismatch_edges.list[0].out_node;
+            **/
+
             /***********需要注释掉********* */
             if (nodeID != currentNodeID + 1)
             {
@@ -1699,6 +1538,11 @@ long long type, long long edgeID)
     else if (type == DELETION)
     {
         nodeID = backbone->g_nodes.list[currentNodeID].deletion_edges.list[edgeID].out_node;
+        dumy->corrected_base += nodeID - currentNodeID;
+        ///currentNodeID = i means self_string[i - 1]
+        add_cigar_record(self_string + currentNodeID, nodeID - currentNodeID, current_cigar, DELETION);
+
+
 
         /***********需要注释掉********* */
         if (!(nodeID >= backbone->s_start_nodeID && nodeID <= backbone->s_end_nodeID))
@@ -1711,7 +1555,7 @@ long long type, long long edgeID)
         }
         
         /***********需要注释掉********* */
-        dumy->corrected_base += nodeID - currentNodeID;
+        
         return nodeID;
     }
     else if (type == INSERTION)
@@ -1725,17 +1569,22 @@ long long type, long long edgeID)
         for (i = 0; i < step; i++)
         {
             add_base_to_correct_read_directly(dumy, backbone->g_nodes.list[nodeID].base);
+            add_cigar_record(&backbone->g_nodes.list[nodeID].base, 1, current_cigar, INSERTION);
             ///只有一条边
             nodeID = backbone->g_nodes.list[nodeID].insertion_edges.list[0].out_node;
         }
         dumy->corrected_base += step;
 
+        ///currentNodeID = i means self_string[i - 1]
+        ///add_cigar_record(self_string + currentNodeID, step, current_cigar, INSERTION);
+
         /***********需要注释掉********* */
         if (nodeID != currentNodeID)
         {
             fprintf(stderr, "error insertion\n");
         }
         /***********需要注释掉********* */
+
         return nodeID;
     }
     else
@@ -1745,7 +1594,7 @@ long long type, long long edgeID)
     
 }
 
-void get_seq_from_Graph(Graph* backbone, Correct_dumy* dumy)
+void get_seq_from_Graph(Graph* backbone, Correct_dumy* dumy, Cigar_record* current_cigar, char* self_string)
 {
     long long new_seq_length = 0;
     long long currentNodeID;
@@ -1839,12 +1688,15 @@ void get_seq_from_Graph(Graph* backbone, Correct_dumy* dumy)
             ///这种情况下矫正
             if(max_count >= total_count*CORRECT_THRESHOLD)
             {
-                currentNodeID = add_path_to_correct_read(backbone, dumy, currentNodeID, max_type, max_edge);
+                currentNodeID = add_path_to_correct_read(backbone, dumy, currentNodeID, max_type, max_edge, current_cigar,
+                self_string);
             }
             else  ///不矫正, 直接取下一个backbone节点
             {
                 currentNodeID++;
                 add_base_to_correct_read_directly(dumy, backbone->g_nodes.list[currentNodeID].base);
+
+                add_cigar_record(&(backbone->g_nodes.list[currentNodeID].base), 1, current_cigar, 0);
             }
             
 
@@ -2296,7 +2148,7 @@ void get_seq_from_Graph_Len2(Graph* backbone, long long backbone_start, long lon
 
 
 void window_consensus(char* r_string, long long window_start, long long window_end, 
-overlap_region_alloc* overlap_list, Correct_dumy* dumy, All_reads* R_INF, Graph* g)
+overlap_region_alloc* overlap_list, Correct_dumy* dumy, All_reads* R_INF, Graph* g, Cigar_record* current_cigar)
 {
     clear_Graph(g);
 
@@ -2358,7 +2210,7 @@ overlap_region_alloc* overlap_list, Correct_dumy* dumy, All_reads* R_INF, Graph*
                     y_string, y_length, &(overlap_list->list[overlapID].w_list[windowID].cigar), startNodeID, endNodeID);
     }
 
-    get_seq_from_Graph(g, dumy);
+    get_seq_from_Graph(g, dumy, current_cigar, backbone);
 
 
     ///get_seq_from_Graph(g, startNodeID, endNodeID, dumy);
@@ -2444,9 +2296,9 @@ overlap_region_alloc* overlap_list, Correct_dumy* dumy, All_reads* R_INF, Graph*
 
 
 void generate_consensus(overlap_region_alloc* overlap_list, All_reads* R_INF, 
-                        UC_Read* g_read, Correct_dumy* dumy, Graph* g)
+                        UC_Read* g_read, Correct_dumy* dumy, Graph* g, Cigar_record* current_cigar)
 {
-
+    clear_Cigar_record(current_cigar);
     long long window_num = (g_read->length + WINDOW - 1) / WINDOW;
     long long i, j, overlap_length;
     long long window_start, window_end;
@@ -2468,6 +2320,7 @@ void generate_consensus(overlap_region_alloc* overlap_list, All_reads* R_INF,
         ///flag返回的是重叠数量
         ///dumy->length返回的是有效完全重叠的数量
         ///dumy->lengthNT返回的是有效不完全重叠的数量
+        ///return overlaps that is overlaped with [window_start, window_end]
         flag = get_available_interval(window_start, window_end, overlap_list, dumy);
         switch (flag)
         {
@@ -2488,11 +2341,12 @@ void generate_consensus(overlap_region_alloc* overlap_list, All_reads* R_INF,
         ///重叠窗口数，也就是coverage大小
         if(dumy->length >= MIN_COVERAGE_THRESHOLD)
         {
-            window_consensus(g_read->seq, window_start, window_end, overlap_list, dumy, R_INF, g);
+            window_consensus(g_read->seq, window_start, window_end, overlap_list, dumy, R_INF, g, current_cigar);
         }
         else
         {
             add_segment_to_correct_read(dumy, g_read->seq + window_start, window_end - window_start + 1);
+            add_cigar_record(g_read->seq + window_start, window_end - window_start + 1, current_cigar, 0);
         }
         
         window_start = window_start + WINDOW;
@@ -2507,6 +2361,7 @@ void generate_consensus(overlap_region_alloc* overlap_list, All_reads* R_INF,
     if (window_start < g_read->length)
     {
         add_segment_to_correct_read(dumy, g_read->seq + window_start, g_read->length - window_start);
+        add_cigar_record(g_read->seq + window_start, g_read->length - window_start, current_cigar, 0);
     }
     
 
@@ -2533,11 +2388,1457 @@ void generate_consensus(overlap_region_alloc* overlap_list, All_reads* R_INF,
     // }
     /***********************要注释掉*************************/
 }
+
+
+
+
+
+
+
+
+
+void markSNP(
+long long window_offset,
+long long x_total_start, long long x_length, 
+long long y_total_start, long long y_length, 
+CIGAR* cigar, haplotype_evdience_alloc* hap)
+{
+    
+    int x_i, y_i, cigar_i;
+    x_i = 0;
+    y_i = 0;
+    cigar_i = 0;
+    int operation;
+    int operationLen;
+    int i;
+    long long inner_offset = x_total_start - window_offset;
+
+    
+    ///note that node 0 is the start node
+    ///0 is match, 1 is mismatch, 2 is up, 3 is left
+    ///2是x缺字符（y多字符），而3是y缺字符（x多字符）
+    while (cigar_i < cigar->length)
+    {
+        operation = cigar->C_C[cigar_i];
+        operationLen = cigar->C_L[cigar_i];
+
+        ///这种情况代表匹配和mismatch
+        if (operation == 0)
+        {
+            x_i += operationLen;
+            y_i += operationLen;
+        }
+        else if(operation == 1)
+        {
+            for (i = 0; i < operationLen; i++)
+            {
+                /**
+                if(inner_offset + x_i >= WINDOW)
+                {
+                    fprintf(stderr, "error\n");
+                }
+                **/
+                hap->flag[inner_offset + x_i]++;
+                x_i++;
+                y_i++;
+            }
+        }///insertion
+        else if (operation == 2)
+        {
+            y_i += operationLen;
+        }
+        else if (operation == 3)
+        {
+            x_i += operationLen;
+        }
+        
+        cigar_i++;
+    }
+    /**
+    if(x_i != x_length || y_i != y_length)
+    {
+        fprintf(stderr, "x_i: %d, x_length: %d\n", x_i, x_length);
+        fprintf(stderr, "y_i: %d, y_length: %d\n", y_i, y_length);
+    }
+    
+   int no_zero = 0;
+   for (i = 0; i < WINDOW; i++)
+   {
+       if(hap->flag[i] != 0)
+       {
+           no_zero++;
+       }
+   }
+
+   fprintf(stderr, "no_zero: %d\n", no_zero);
+   **/
+   
+}
+
+
+
+
+void addSNPtohaplotype(
+long long window_offset, int overlapID,
+char* x_string, long long x_total_start, long long x_length, 
+char* y_string, long long y_total_start, long long y_length, 
+CIGAR* cigar, haplotype_evdience_alloc* hap)
+{
+    
+    int x_i, y_i, cigar_i;
+    x_i = 0;
+    y_i = 0;
+    cigar_i = 0;
+    int operation;
+    int operationLen;
+    int i;
+    long long inner_offset = x_total_start - window_offset;
+    haplotype_evdience ev;
+    
+    ///note that node 0 is the start node
+    ///0 is match, 1 is mismatch, 2 is up, 3 is left
+    ///2是x缺字符（y多字符），而3是y缺字符（x多字符）
+    while (cigar_i < cigar->length)
+    {
+        operation = cigar->C_C[cigar_i];
+        operationLen = cigar->C_L[cigar_i];
+
+        ///这种情况代表匹配和mismatch
+        if (operation == 0)
+        {
+            for (i = 0; i < operationLen; i++)
+            {
+                if(hap->flag[inner_offset] > FLAG_THRE)
+                {
+                    ev.misBase =  y_string[y_i];
+                    ev.overlapID = overlapID;
+                    ev.site = x_total_start + x_i;
+                    ev.overlapSite = y_total_start + y_i;
+                    ev.type = 0;
+                    addHaplotypeEvdience(hap, &ev);
+                }
+
+
+                inner_offset++;
+                x_i++;
+                y_i++;
+            }
+
+        }
+        else if(operation == 1)
+        {
+            for (i = 0; i < operationLen; i++)
+            {
+
+                if(hap->flag[inner_offset] > FLAG_THRE)
+                {
+                    ev.misBase =  y_string[y_i];
+                    ev.overlapID = overlapID;
+                    ev.site = x_total_start + x_i;
+                    ev.overlapSite = y_total_start + y_i;
+                    ev.type = 1;
+                    addHaplotypeEvdience(hap, &ev);
+                }
+
+                inner_offset++;
+                x_i++;
+                y_i++;
+            }
+        }///insertion
+        else if (operation == 2)
+        {
+            y_i += operationLen;
+        }
+        else if (operation == 3)
+        {
+
+            for (i = 0; i < operationLen; i++)
+            {
+                if(hap->flag[inner_offset] > FLAG_THRE)
+                {
+                    ev.misBase =  'N';
+                    ev.overlapID = overlapID;
+                    ev.site = x_total_start + x_i;
+                    ev.overlapSite = y_total_start + y_i;
+                    ev.type = 2;
+                    addHaplotypeEvdience(hap, &ev);
+                }
+
+                inner_offset++;
+                x_i++;
+            }
+        }
+        
+        cigar_i++;
+    }
+}
+
+
+
+
+
+
+
+
+
+
+
+
+void cluster(char* r_string, long long window_start, long long window_end, 
+overlap_region_alloc* overlap_list, Correct_dumy* dumy, All_reads* R_INF, haplotype_evdience_alloc* hap)
+{
+    long long x_start;
+    long long x_length; 
+    char* x_string;
+    char* y_string;
+    long long i;
+    long long y_start, y_length;
+    long long overlapID, windowID;
+    long long startNodeID, endNodeID, currentNodeID;
+
+    RsetInitHaplotypeEvdienceFlag(hap);
+
+    long long correct_x_pos_s;
+    long long inner_window_offset;
+
+
+    ///与当前window重叠的所有overlap
+    ///first mark all snp pos
+    for (i = 0; i < dumy->length; i++)
+    {
+        ///这个是那个overlap的ID，而不是overlap里对应窗口的ID
+        overlapID = dumy->overlapID[i];
+
+        ///overlap_list->list[overlapID].x_pos_s is the begining of the whole overlap
+        correct_x_pos_s = (overlap_list->list[overlapID].x_pos_s / WINDOW) * WINDOW;
+        ///window_start is the begining of this window in the whole x_read
+        windowID = (window_start - correct_x_pos_s) / WINDOW;
+
+        ///如果这个window不匹配，跳过
+        if (overlap_list->list[overlapID].w_list[windowID].y_end == -1)
+        {
+            continue;
+        }
+        
+        ///both x_start and y_start are the begining of the whole x_read and y_read
+        x_start = overlap_list->list[overlapID].w_list[windowID].x_start;
+        x_length = overlap_list->list[overlapID].w_list[windowID].x_end 
+                - overlap_list->list[overlapID].w_list[windowID].x_start + 1;
+
+        y_start = overlap_list->list[overlapID].w_list[windowID].y_start;
+        y_length = overlap_list->list[overlapID].w_list[windowID].y_end
+                - overlap_list->list[overlapID].w_list[windowID].y_start + 1;
+
+
+        markSNP(window_start, x_start, x_length, y_start, y_length, &(overlap_list->list[overlapID].w_list[windowID].cigar), 
+        hap);
+    }
+
+
+    for (i = 0; i < WINDOW; i++)
+    {
+        if(hap->flag[i] > FLAG_THRE)
+        {
+            hap->snp++;
+        }
+    }
+
+
+
+    ///add the information related to snp to haplotype_evdience_alloc
+    for (i = 0; i < dumy->length; i++)
+    {
+        ///这个是那个overlap的ID，而不是overlap里对应窗口的ID
+        overlapID = dumy->overlapID[i];
+
+        ///overlap_list->list[overlapID].x_pos_s is the begining of the whole overlap
+        correct_x_pos_s = (overlap_list->list[overlapID].x_pos_s / WINDOW) * WINDOW;
+        ///window_start is the begining of this window in the whole x_read
+        windowID = (window_start - correct_x_pos_s) / WINDOW;
+
+        ///如果这个window不匹配，跳过
+        if (overlap_list->list[overlapID].w_list[windowID].y_end == -1)
+        {
+            continue;
+        }
+        
+        ///both x_start and y_start are the begining of the whole x_read and y_read
+        x_start = overlap_list->list[overlapID].w_list[windowID].x_start;
+        x_length = overlap_list->list[overlapID].w_list[windowID].x_end 
+                - overlap_list->list[overlapID].w_list[windowID].x_start + 1;
+
+        y_start = overlap_list->list[overlapID].w_list[windowID].y_start;
+        y_length = overlap_list->list[overlapID].w_list[windowID].y_end
+                - overlap_list->list[overlapID].w_list[windowID].y_start + 1;
+
+
+        recover_UC_Read_sub_region(dumy->overlap_region, y_start, y_length, overlap_list->list[overlapID].y_pos_strand, 
+                R_INF, overlap_list->list[overlapID].y_id);
+
+        x_string = r_string + x_start;
+        y_string = dumy->overlap_region;
+
+
+        addSNPtohaplotype(window_start, overlapID, x_string, x_start, x_length, 
+        y_string, y_start, y_length, &(overlap_list->list[overlapID].w_list[windowID].cigar), 
+        hap);        
+    }
+
+
+
+}
+
+int cmp_haplotype_evdience(const void * a, const void * b)
+{
+    if ((*(haplotype_evdience*)a).site != (*(haplotype_evdience*)b).site)
+    {
+        return (*(haplotype_evdience*)a).site > (*(haplotype_evdience*)b).site ? 1 : -1; 
+    }
+    else
+    {
+        if ((*(haplotype_evdience*)a).type != (*(haplotype_evdience*)b).type)
+        {
+            return (*(haplotype_evdience*)a).type > (*(haplotype_evdience*)b).type ? 1 : -1; 
+        }
+        else
+        {
+            if ((*(haplotype_evdience*)a).misBase != (*(haplotype_evdience*)b).misBase)
+            {
+                return (*(haplotype_evdience*)a).misBase > (*(haplotype_evdience*)b).misBase ? 1 : -1; 
+            }
+            else
+            {
+                return 0;
+            }
+            
+        }
+        
+    }
+    
+    
+}
+
+
+int cmp_snp_stats(const void * a, const void * b)
+{
+    if ((*(SnpStats*)a).score != (*(SnpStats*)b).score)
+    {
+        return (*(SnpStats*)a).score < (*(SnpStats*)b).score ? 1 : -1; 
+    }
+    else
+    {
+        if ((*(SnpStats*)a).occ_2 != (*(SnpStats*)b).occ_2)
+        {
+            return (*(SnpStats*)a).occ_2 > (*(SnpStats*)b).occ_2 ? 1 : -1; 
+        }
+        else
+        {
+            return 0;
+        }   
+    }
+}
+
+void debug_hap_information(overlap_region_alloc* overlap_list, All_reads* R_INF, 
+                        UC_Read* g_read, haplotype_evdience_alloc* hap,
+                        Correct_dumy* dumy)
+{
+    int i, overlapID, y_ID, y_Strand;
+    long long x_start;
+    long long x_length; 
+    char* x_string;
+    char* y_string;
+    long long y_start;
+    long long y_length;
+
+    for (i = 0; i < hap->length; i++)
+    {
+        if(hap->list[i].type < 2)
+        {
+            overlapID = hap->list[i].overlapID;
+            x_start = hap->list[i].site; 
+            y_start = hap->list[i].overlapSite;
+
+
+            y_ID = overlap_list->list[overlapID].y_id;
+            y_Strand = overlap_list->list[overlapID].y_pos_strand;
+
+            recover_UC_Read_sub_region(dumy->overlap_region, y_start, 1, y_Strand, R_INF, y_ID);
+
+            x_string = g_read->seq + x_start;
+            y_string = dumy->overlap_region;
+
+            if(y_string[0] != hap->list[i].misBase)
+            {
+                fprintf(stderr, "y_string[0]: %c, hap->list[i].misBase: %c\n",
+                    y_string[0], hap->list[i].misBase);
+            }
+
+
+            if(hap->list[i].type == 0)
+            {
+                if(x_string[0] != y_string[0])
+                {
+                    fprintf(stderr, "x_string[0]: %c, y_string[0]: %c\n",
+                    x_string[0], y_string[0]);
+
+
+                }
+            }
+            else if(hap->list[i].type == 0)
+            {
+                if(x_string[0] == y_string[0])
+                {
+                    fprintf(stderr, "x_string[0]: %c, y_string[0]: %c\n",
+                    x_string[0], y_string[0]);
+                }
+
+            }
+        }
+
+        
+
+    }
+}
+
+
+
+
+
+
+
+
+
+int debug_split_sub_list(haplotype_evdience_alloc* hap, 
+haplotype_evdience* sub_list, long long sub_length, long long num_haplotype)
+{
+    long long i = 0;
+    long long occ_0 = 0;
+    long long occ_1 = 0;
+    long long occ_1_array[5];
+    memset(occ_1_array, 0, sizeof(long long) * 5);
+    long long occ_2 = 0;
+
+
+    for (i = 0; i < sub_length; i++)
+    {
+        if(sub_list[i].type == 0)
+        {
+            occ_0++;
+        }
+        else if(sub_list[i].type == 1)
+        {
+            occ_1_array[seq_nt6_table[(uint8_t)(sub_list[i].misBase)]]++;
+            occ_1++;
+        }
+        else if(sub_list[i].type == 2)
+        {
+            occ_2++;
+        }
+    }
+
+    /**
+     1. if occ_0 = 0, that means all overlaps are different with this read at this site
+     2. it is not possible that occ_1 = 0,
+     3. if occ_1 = 1, there are only one difference. It must be a sequencing error.
+    **/
+    if(occ_0 == 0 || occ_1 <= 1)
+    {
+        return 0;
+    }
+
+    ///note: if the max value except type0 is type2
+    ///that means this is no snp hapolyte
+    long long max = occ_2;
+    long long max_i = -1;
+
+    for (i = 0; i < 5; i++)
+    {
+        if(occ_1_array[i] > max)
+        {
+            max = occ_1_array[i];
+            max_i = i;
+        }
+    }
+
+
+
+    if(max_i == -1)
+    {
+        return 0;
+    }
+
+    if(max <= 1)
+    {
+        return 0;
+    }
+
+    ///if we have two max
+    for (i = 0; i < 5; i++)
+    {
+        if(occ_1_array[i] == max && i != max_i)
+        {
+            return 0;
+        }
+    }
+
+    long long new_0 = occ_0 + 1;
+    long long new_total = sub_length + 1;
+    ///note: here occ_0++ since the read itself has a type0
+    double available = new_0 + max;
+    double threshold = 0.95;
+    available = available/((double)(new_total));
+    if(available < threshold)
+    {
+        return 0;
+    }
+
+    ///if we just have one snp, we need to phase it carefully
+    if(num_haplotype == 1)
+    {
+        ///we must have just 1 match and 1 mismatch
+        ///any other types are not good
+        if(new_0 + max != new_total)
+        {
+            return 0;
+        }
+        
+        if(filter_snp(new_0, max, new_total) == 0)
+        {
+            return 0;
+        }
+    }
+    /**
+    if(filter_snp(new_0, max, new_total) == 0)
+    {
+        return 0;
+    }
+    **/
+
+    
+
+
+    ///for each calculated snp, find if it is at snp matrix
+    for (i = 0; i < hap->available_snp; i++)
+    {
+        if(hap->snp_stat[i].site == sub_list[0].site)
+        {
+
+            
+            int j = 0;
+            int vectorID = hap->snp_stat[i].id;
+            int8_t* vector = Get_SNP_Vector((*hap), vectorID);
+
+            if(hap->snp_stat[i].occ_0 != occ_0)
+            {
+                fprintf(stderr, "error occ0\n");
+            }
+
+            if(hap->snp_stat[i].occ_1 != occ_1_array[max_i])
+            {
+                fprintf(stderr, "error occ1\n");
+            }
+
+            if(hap->snp_stat[i].overlap_num != sub_length)
+            {
+                fprintf(stderr, "error overlap_num\n");
+            }
+
+            if(hap->snp_stat[i].overlap_num != hap->snp_stat[i].occ_0 + 
+            hap->snp_stat[i].occ_1 + hap->snp_stat[i].occ_2)
+            {
+                fprintf(stderr, "error overlap_num\n");
+            }
+
+            ///for each element in snp vector, find if it is in calculated dataset
+            for (j = 0; j < Get_SNP_Vector_Length((*hap)); j++)
+            {
+                if(vector[j] != -1)
+                {
+                    int x_i = 0;
+                    for (x_i = 0; x_i < sub_length; x_i++)
+                    {
+                        if(j == sub_list[x_i].overlapID)
+                        {
+                            break;
+                        }
+                    }
+
+                    if(x_i == sub_length)
+                    {
+                        fprintf(stderr, "error: j: %d\n",j);
+                    }
+                    else
+                    {
+                        if(vector[j] == 0 || sub_list[x_i].type == 0)
+                        {
+                            if(vector[j] != sub_list[x_i].type)
+                            {
+                                fprintf(stderr, "error: 0: %d\n",j);
+                            }
+                        }
+
+                        if(vector[j] == 1)
+                        {
+                            if(sub_list[x_i].type != 1)
+                            {
+                                fprintf(stderr, "-error: 1: %d\n",j);
+                            }
+
+
+                            if(sub_list[x_i].type == 1 && sub_list[x_i].misBase != s_H[max_i])
+                            {
+                                fprintf(stderr, "+error: 1: %d\n",j);
+                            }
+                        }
+
+                        if(vector[j] == 2)
+                        {
+                            if(sub_list[x_i].type != 2)
+                            {
+                                if(sub_list[x_i].type == 1 && sub_list[x_i].misBase != s_H[max_i])
+                                {
+                                    ;
+                                }
+                                else
+                                {
+                                    fprintf(stderr, "error: 2: %d\n",j);
+                                }
+                    
+                            }
+
+                        }
+                        
+         
+                    }
+                    
+                }
+            }
+
+
+            ///for each calculated data, find if it is in snp vector
+            for (j = 0; j < sub_length; j++)
+            {
+                if(vector[sub_list[j].overlapID] != sub_list[j].type)
+                {
+
+                    if(vector[sub_list[j].overlapID] == 2 && sub_list[j].type == 1 && sub_list[j].misBase != s_H[max_i])
+                    {
+                        ;
+                    }
+                    else
+                    {
+                            fprintf(stderr, "vector[sub_list[j].site]: %d, sub_list[j].type: %d\n",
+                        vector[sub_list[j].overlapID], sub_list[j].type);
+                    }
+                }
+            }
+            
+            
+            break;
+        }
+    }
+
+    
+    if(i == hap->available_snp)
+    {
+        fprintf(stderr, "error\n");
+    }    
+    
+
+    /**
+    fprintf(stderr, "new_0: %d, occ_0: %d, max: %d, max_i: %d, sub_length: %d, new_total: %d, available: %lf\n", 
+    new_0, occ_0, max, max_i, sub_length, new_total, available);
+    for (i = 0; i < sub_length; i++)
+    {
+        
+        fprintf(stderr, "i: %d, site: %d, type: %d, char: %c, ID: %d, name: %.*s\n", 
+    i, sub_list[i].site, sub_list[i].type, sub_list[i].misBase, sub_list[i].overlapID, 
+    Get_NAME_LENGTH((*R_INF), overlap_list->list[sub_list[i].overlapID].y_id), 
+    Get_NAME((*R_INF),overlap_list->list[sub_list[i].overlapID].y_id)); 
+
+    }
+    fprintf(stderr, "\n");
+    **/
+    
+    
+    return 1;
+    
+
+    
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+int debug_snp_matrix(haplotype_evdience_alloc* hap)
+{
+    uint64_t pre_site = (uint64_t)-1;
+    uint64_t num_of_snps = 0;
+    long long pre_i = -1;
+    long long sub_length;
+    haplotype_evdience* sub_list;
+    long long i;
+    long long a_snp = 0;
+
+
+    ////split reads
+    for (i = 0; i < hap->length; i++)
+    {
+        if(pre_site != hap->list[i].site)
+        {
+            if(i != 0)
+            {
+                sub_list = hap->list + pre_i;
+                sub_length = i - pre_i;
+                ///debug_total_length = debug_total_length + sub_length;
+                a_snp += debug_split_sub_list(hap, sub_list, sub_length, hap->snp);
+            }
+            num_of_snps++;
+            pre_site = hap->list[i].site;
+            pre_i = i;
+        }
+    }
+
+    if(pre_i != -1)
+    {
+        sub_list = hap->list + pre_i;
+        sub_length = i - pre_i;
+        ///debug_total_length = debug_total_length + sub_length;
+        a_snp += debug_split_sub_list(hap, sub_list, sub_length, hap->snp);
+    }
+
+    
+    if(a_snp != hap->available_snp)
+    {
+        fprintf(stderr, "a_snp: %d, available_snp: %d\n", 
+        a_snp, hap->available_snp);
+    }
+    
+
+}
+
+int split_sub_list(haplotype_evdience_alloc* hap, 
+haplotype_evdience* sub_list, long long sub_length, long long num_haplotype, 
+overlap_region_alloc* overlap_list, All_reads* R_INF)
+{
+    long long i = 0;
+    long long occ_0 = 0;
+    long long occ_1 = 0;
+    long long occ_1_array[5];
+    memset(occ_1_array, 0, sizeof(long long) * 5);
+    long long occ_2 = 0;
+
+
+    for (i = 0; i < sub_length; i++)
+    {
+        if(sub_list[i].type == 0)
+        {
+            occ_0++;
+        }
+        else if(sub_list[i].type == 1)
+        {
+            occ_1_array[seq_nt6_table[(uint8_t)(sub_list[i].misBase)]]++;
+            occ_1++;
+        }
+        else if(sub_list[i].type == 2)
+        {
+            occ_2++;
+        }
+    }
+
+    /**
+     1. if occ_0 = 0, that means all overlaps are different with this read at this site
+     2. it is not possible that occ_1 = 0,
+     3. if occ_1 = 1, there are only one difference. It must be a sequencing error.
+    **/
+    if(occ_0 == 0 || occ_1 <= 1)
+    {
+        return 0;
+    }
+
+    ///note: if the max value except type0 is type2
+    ///that means this is no snp hapolyte
+    long long max = occ_2;
+    long long max_i = -1;
+
+    for (i = 0; i < 5; i++)
+    {
+        if(occ_1_array[i] > max)
+        {
+            max = occ_1_array[i];
+            max_i = i;
+        }
+    }
+
+
+
+    if(max_i == -1)
+    {
+        return 0;
+    }
+
+    if(max <= 1)
+    {
+        return 0;
+    }
+
+    ///if we have two max
+    for (i = 0; i < 5; i++)
+    {
+        if(occ_1_array[i] == max && i != max_i)
+        {
+            return 0;
+        }
+    }
+
+    long long new_0 = occ_0 + 1;
+    long long new_total = sub_length + 1;
+    ///note: here occ_0++ since the read itself has a type0
+    double available = new_0 + max;
+    double threshold = 0.95;
+    available = available/((double)(new_total));
+    if(available < threshold)
+    {
+        return 0;
+    }
+
+    ///if we just have one snp, we need to phase it carefully
+    if(num_haplotype == 1)
+    {
+        ///we must have just 1 match and 1 mismatch
+        ///any other types are not good
+        if(new_0 + max != new_total)
+        {
+            return 0;
+        }
+
+        if(filter_snp(new_0, max, new_total) == 0)
+        {
+            return 0;
+        }
+    }
+
+    /**
+    if(filter_snp(new_0, max, new_total) == 0)
+    {
+        return 0;
+    }
+    **/
+    
+        
+
+
+    
+
+    
+    InsertSNPVector(hap, sub_list, sub_length, s_H[max_i]);
+    
+    
+
+    /**
+    fprintf(stderr, "new_0: %d, occ_0: %d, max: %d, max_i: %d, sub_length: %d, new_total: %d, available: %lf\n", 
+    new_0, occ_0, max, max_i, sub_length, new_total, available);
+    for (i = 0; i < sub_length; i++)
+    {
+        
+        fprintf(stderr, "i: %d, site: %d, type: %d, char: %c, ID: %d, name: %.*s\n", 
+    i, sub_list[i].site, sub_list[i].type, sub_list[i].misBase, sub_list[i].overlapID, 
+    Get_NAME_LENGTH((*R_INF), overlap_list->list[sub_list[i].overlapID].y_id), 
+    Get_NAME((*R_INF),overlap_list->list[sub_list[i].overlapID].y_id)); 
+
+    }
+    fprintf(stderr, "\n");
+    **/
+    
+    
+    
+    
+
+    
+}
+
+
+int calculate_distance_snp_vector(int8_t *vector1, int8_t *vector2, int Len)
+{
+    int i;
+    for (i = 0; i < Len; i++)
+    {
+        if(vector1[i] != vector2[i])
+        {
+            if ((vector1[i] == 0 || vector1[i] == 1) && (vector2[i] == 0 || vector2[i] == 1))
+            {
+                return 1;
+            }
+        }
+    }
+    
+    return 0;
+}
+
+void print_core_snp(haplotype_evdience_alloc* hap)
+{
+    int i, j;
+    for (i = 0; i < hap->core_snp; i++)
+    {
+        fprintf(stderr, "core(i): %d, site: %d, occ_0: %d, occ_1: %d, occ_2: %d, score: %d\n", 
+    i, hap->snp_stat[i].site, hap->snp_stat[i].occ_0, hap->snp_stat[i].occ_1,
+    hap->snp_stat[i].occ_2,
+    hap->snp_stat[i].score); 
+
+        int vectorID = hap->snp_stat[i].id;
+        int8_t* vector = Get_SNP_Vector((*hap), vectorID);
+
+        for (j = 0; j < hap->overlap; j++)
+        {
+            if(vector[j] == 0)
+            {
+                fprintf(stderr, "type: %d, ID: %d\n", vector[j], j);
+            }
+        }
+
+
+        for (j = 0; j < hap->overlap; j++)
+        {
+            if(vector[j] == 1)
+            {
+                fprintf(stderr, "type: %d, ID: %d\n", vector[j], j);
+            }
+        }
+
+
+        for (j = 0; j < hap->overlap; j++)
+        {
+            if(vector[j] == 2)
+            {
+                fprintf(stderr, "type: %d, ID: %d\n", vector[j], j);
+            }
+        }
+        
+    }
+}
+
+void merge_snp_vectors(haplotype_evdience_alloc* hap, int diff_vector_ID)
+{
+    int8_t *r_vector = Get_Result_SNP_Vector((*hap));
+    int vectorLen = Get_SNP_Vector_Length((*hap));
+    memset(r_vector, -1, vectorLen);
+    hap->result_stat.occ_0 = 0;
+    hap->result_stat.occ_1 = 0;
+
+    int8_t* vector;
+    int vectorID;
+    int i, j;
+
+    for (i = 0; i < hap->core_snp; i++)
+    {
+        if(i == diff_vector_ID)
+        {
+            continue;
+        }
+
+        vectorID = hap->snp_stat[i].id;
+        vector = Get_SNP_Vector((*hap), vectorID);
+
+        for (j = 0; j < vectorLen; j++)
+        {
+            if(r_vector[j] == -1)
+            {
+                if(vector[j] == 0)
+                {
+                    hap->result_stat.occ_0++;
+                    r_vector[j] = vector[j];
+                }
+                else if(vector[j] == 1)
+                {
+                    hap->result_stat.occ_1++;
+                    r_vector[j] = vector[j];
+                }
+            }
+            else ///can debug here
+            {
+                if((vector[j] != -1 && vector[j] != 2 && vector[j] != r_vector[j]))
+                {
+                    fprintf(stderr, "j: %d, vector[j]: %d, r_vector[j]: %d, hap->core_snp: %d, diff_vector_ID: %d\n",
+                    j, vector[j], r_vector[j], hap->core_snp, diff_vector_ID);
+
+                    print_core_snp(hap);
+                }
+            }
+            
+            
+            
+        }
+    }
+
+    hap->result_stat.overlap_num = hap->result_stat.occ_0 + hap->result_stat.occ_1;
+}
+
+void add_to_result_snp_vector(haplotype_evdience_alloc* hap, int8_t *new_vector, int Len)
+{
+    int8_t *r_vector = Get_Result_SNP_Vector((*hap));
+    int j;
+    for (j = 0; j < Len; j++)
+    {
+        if(r_vector[j] == -1)
+        {
+            if(new_vector[j] == 0)
+            {
+                hap->result_stat.occ_0++;
+                r_vector[j] = new_vector[j];
+            }
+            else if(new_vector[j] == 1)
+            {
+                hap->result_stat.occ_1++;
+                r_vector[j] = new_vector[j];
+            }
+        }
+        ///can debug here
+    }
+
+    hap->result_stat.overlap_num = hap->result_stat.occ_0 + hap->result_stat.occ_1;
+}
+
+
+
+
+int merge_snp_vectors_and_test(haplotype_evdience_alloc* hap, int diff_vector_ID)
+{
+    int8_t *r_vector = Get_Result_SNP_Vector((*hap));
+    int vectorLen = Get_SNP_Vector_Length((*hap));
+    memset(r_vector, -1, vectorLen);
+    hap->result_stat.occ_0 = 0;
+    hap->result_stat.occ_1 = 0;
+
+    int8_t* vector;
+    int vectorID;
+    int i, j;
+
+    for (i = 0; i < hap->core_snp; i++)
+    {
+        if(i == diff_vector_ID)
+        {
+            continue;
+        }
+
+        vectorID = hap->snp_stat[i].id;
+        vector = Get_SNP_Vector((*hap), vectorID);
+
+        for (j = 0; j < vectorLen; j++)
+        {
+            if(r_vector[j] == -1)
+            {
+                if(vector[j] == 0)
+                {
+                    hap->result_stat.occ_0++;
+                    r_vector[j] = vector[j];
+                }
+                else if(vector[j] == 1)
+                {
+                    hap->result_stat.occ_1++;
+                    r_vector[j] = vector[j];
+                }
+            }
+            else ///can debug here
+            {
+                ///has confilict
+                if(vector[j] != -1 && vector[j] != 2 && vector[j] != r_vector[j])
+                {
+                    return 0;
+                }
+            }
+        }
+    }
+
+    hap->result_stat.overlap_num = hap->result_stat.occ_0 + hap->result_stat.occ_1;
+
+    return 1;
+}
+
+int generate_haplotypes(haplotype_evdience_alloc* hap)
+{
+    int j;
+
+    int vectorID, vectorID2;
+    int diff_core_vector = 0;
+    int diff_vector_ID = -1;
+    int8_t *vector, *vector2;
+    
+    if(hap->core_snp == 0)
+    {
+        return 0;
+    }
+
+    ///the hap->core_snp is used to find centriod
+
+            
+    ///if there are <5 vectors in core_snp, we didn't allow different vector
+    if (hap->core_snp < 5)
+    {
+        if(merge_snp_vectors_and_test(hap, -1) == 0)
+        {
+            return 0;
+        }
+    }
+    else ///for vectors in core_snp, we allow at most one different vector when there are >= 5 vectors in core_snp
+    {
+        ///there are two condition: 1. vector 0 is the different one. 2. vector 0 is not the different one
+
+        
+        diff_vector_ID = -1;
+        ///first try to merge all vector together
+        if(merge_snp_vectors_and_test(hap, -1) == 0)
+        {
+            for (j = hap->core_snp - 1; j >= 0; j--)
+            {
+                if(merge_snp_vectors_and_test(hap, j) == 1)
+                {
+                    diff_vector_ID = j;
+                    break;
+                }
+            }
+
+            if(j == -1)
+            {
+                return 0;
+            }
+        }
+    }
+
+
+    
+
+    ///after merge, we get result vector
+    vector = Get_Result_SNP_Vector((*hap));
+    ///and for each non-core snp vector, if it has no conflict with result vector
+    /// add it to result vector
+    for (j = hap->core_snp; j < hap->available_snp; j++)
+    {
+        vectorID2 = hap->snp_stat[j].id;
+        vector2 = Get_SNP_Vector((*hap), vectorID2);
+        if(calculate_distance_snp_vector(vector, vector2, Get_SNP_Vector_Length((*hap))) == 0)
+        {
+            add_to_result_snp_vector(hap, vector2, Get_SNP_Vector_Length((*hap)));
+        }
+    }
+
+    ///merge_snp_vectors(hap, diff_vector_ID);
+
+    ///for read only have 1 snp, we need a more strict condition
+    if (hap->core_snp == 1 && 
+    filter_one_snp(hap->result_stat.occ_0 + 1, hap->result_stat.occ_1, 
+    hap->result_stat.overlap_num + 1) == 0)
+    {
+        return 0;
+    }
+
+
+    return 1;
+    
+}
+
+void print_Haplotype(haplotype_evdience_alloc* hap, overlap_region_alloc* overlap_list, All_reads* R_INF)
+{
+        int j, i;
+
+        
+    fprintf(stderr, "\nhap->snp: %d, hap->length: %d, perc: %d, x_name: %.*s\n", 
+    hap->snp, hap->length, (hap->snp == 0? 0: hap->length/hap->snp), 
+    Get_NAME_LENGTH((*R_INF), overlap_list->list[0].x_id), Get_NAME((*R_INF),overlap_list->list[0].x_id)
+    );
+    fprintf(stderr, "hap->available_snp: %d, hap->core_snp:%d\n", 
+    hap->available_snp, hap->core_snp);
+
+    int Len_x, matched_overlap = 0;
+        for (j = 0; j < overlap_list->length; j++)
+        {
+            Len_x = overlap_list->list[j].x_pos_e -  overlap_list->list[j].x_pos_s + 1;
+
+            if (Len_x * OVERLAP_THRESHOLD <=  overlap_list->list[j].align_length)
+            {
+                matched_overlap++;
+            }
+        }
+
+        fprintf(stderr, "occ_0: %d, occ_1: %d, overlap_num: %d, matched_overlap: %d, overlap_list->length: %d\n", 
+    hap->result_stat.occ_0, hap->result_stat.occ_1,
+    hap->result_stat.overlap_num, matched_overlap, overlap_list->length); 
+    
+
+    return;
+
+        fprintf(stderr, "Phaseing sucessfully!\n");
+        fprintf(stderr, "occ_0: %d, occ_1: %d, overlap_num: %d\n", 
+    hap->result_stat.occ_0, hap->result_stat.occ_1,
+    hap->result_stat.overlap_num); 
+
+
+        int8_t* vector = Get_Result_SNP_Vector((*hap));
+        for (j = 0; j < hap->overlap; j++)
+        {
+            if(vector[j] == 0)
+            {
+                fprintf(stderr, "Ptype: %d, ID: %d\n", vector[j], j);
+            }
+        }
+
+
+        for (j = 0; j < hap->overlap; j++)
+        {
+            if(vector[j] == 1)
+            {
+                fprintf(stderr, "Ptype: %d, ID: %d\n", vector[j], j);
+            }
+        }
+
+
+        for (j = 0; j < hap->overlap; j++)
+        {
+            if(vector[j] == 2)
+            {
+                fprintf(stderr, "Ptype: %d, ID: %d\n", vector[j], j);
+            }
+        }
+
+
+    for (i = 0; i < hap->core_snp; i++)
+    {
+        fprintf(stderr, "core(i): %d, site: %d, occ_0: %d, occ_1: %d, occ_2: %d, score: %d\n", 
+    i, hap->snp_stat[i].site, hap->snp_stat[i].occ_0, hap->snp_stat[i].occ_1,
+    hap->snp_stat[i].occ_2,
+    hap->snp_stat[i].score); 
+
+        int vectorID = hap->snp_stat[i].id;
+        int8_t* vector = Get_SNP_Vector((*hap), vectorID);
+
+        for (j = 0; j < hap->overlap; j++)
+        {
+            if(vector[j] == 0)
+            {
+                fprintf(stderr, "type: %d, ID: %d\n", vector[j], j);
+            }
+        }
+
+
+        for (j = 0; j < hap->overlap; j++)
+        {
+            if(vector[j] == 1)
+            {
+                fprintf(stderr, "type: %d, ID: %d\n", vector[j], j);
+            }
+        }
+
+
+        for (j = 0; j < hap->overlap; j++)
+        {
+            if(vector[j] == 2)
+            {
+                fprintf(stderr, "type: %d, ID: %d\n", vector[j], j);
+            }
+        }
+        
+    }
+
+
+
+        for (; i < hap->available_snp; i++)
+    {
+        fprintf(stderr, "i: %d, site: %d, occ_0: %d, occ_1: %d, occ_2: %d, score: %d\n", 
+    i, hap->snp_stat[i].site, hap->snp_stat[i].occ_0, hap->snp_stat[i].occ_1,
+    hap->snp_stat[i].occ_2,
+    hap->snp_stat[i].score); 
+
+        int vectorID = hap->snp_stat[i].id;
+        int8_t* vector = Get_SNP_Vector((*hap), vectorID);
+
+        for (j = 0; j < hap->overlap; j++)
+        {
+            if(vector[j] == 0)
+            {
+                fprintf(stderr, "type: %d, ID: %d\n", vector[j], j);
+            }
+        }
+
+
+        for (j = 0; j < hap->overlap; j++)
+        {
+            if(vector[j] == 1)
+            {
+                fprintf(stderr, "type: %d, ID: %d\n", vector[j], j);
+            }
+        }
+
+
+        for (j = 0; j < hap->overlap; j++)
+        {
+            if(vector[j] == 2)
+            {
+                fprintf(stderr, "type: %d, ID: %d\n", vector[j], j);
+            }
+        }
+        
+    }
+    
+}
+
+void partition_overlaps(overlap_region_alloc* overlap_list, All_reads* R_INF, 
+                        UC_Read* g_read, Correct_dumy* dumy, haplotype_evdience_alloc* hap)
+{
+    ResizeInitHaplotypeEvdience(hap);
+
+    long long window_num = (g_read->length + WINDOW - 1) / WINDOW;
+    long long i, j, overlap_length;
+    long long window_start, window_end;
+
+    long long num_availiable_win = 0;
+    
+
+    window_start = 0;
+    window_end = WINDOW - 1;
+    if (window_end >= g_read->length)
+    {
+        window_end = g_read->length - 1;
+    }
+    int flag;
+    for (i = 0; i < window_num; i++)
+    {
+        dumy->length = 0;
+        dumy->lengthNT = 0;
+        ///flag返回的是重叠数量
+        ///dumy->length返回的是有效完全重叠的数量
+        ///dumy->lengthNT返回的是有效不完全重叠的数量
+        ///return overlaps that is overlaped with [window_start, window_end]
+        flag = get_available_interval(window_start, window_end, overlap_list, dumy);
+        switch (flag)
+        {
+            case 1:    ///找到匹配
+                break;
+            case 0:    ///没找到匹配
+                break;
+            case -2: ///下一个window也不会存在匹配, 直接跳出
+                i = window_num;
+                break;
+        }
+
+
+        ///这个是available overlap里所有window的数量...
+        ///num_availiable_win = num_availiable_win + dumy->length + dumy->lengthNT;
+        num_availiable_win = num_availiable_win + dumy->length;
+
+
+        cluster(g_read->seq, window_start, window_end, overlap_list, dumy, R_INF, hap);
+
+        
+        window_start = window_start + WINDOW;
+        window_end = window_end + WINDOW;
+        if (window_end >= g_read->length)
+        {
+            window_end = g_read->length - 1;
+        }
+        
+    }
+
+
+    
+    ///very time-consuming
+    qsort(hap->list, hap->length, sizeof(haplotype_evdience), cmp_haplotype_evdience);
+
+    ///debug_hap_information(overlap_list, R_INF, g_read, hap, dumy);
+
+    SetSnpMatrix(hap, hap->snp, overlap_list->length);
+
+
+    uint64_t pre_site = (uint64_t)-1;
+    uint64_t num_of_snps = 0;
+    long long pre_i = -1;
+    long long sub_length;
+    haplotype_evdience* sub_list;
+
+    ///long long debug_total_length = 0;
+
+    ////split reads
+    for (i = 0; i < hap->length; i++)
+    {
+        if(pre_site != hap->list[i].site)
+        {
+            if(i != 0)
+            {
+                sub_list = hap->list + pre_i;
+                sub_length = i - pre_i;
+                ///debug_total_length = debug_total_length + sub_length;
+                split_sub_list(hap, sub_list, sub_length, hap->snp, overlap_list, R_INF);
+            }
+            num_of_snps++;
+            pre_site = hap->list[i].site;
+            pre_i = i;
+        }
+    }
+
+    if(pre_i != -1)
+    {
+        sub_list = hap->list + pre_i;
+        sub_length = i - pre_i;
+        ///debug_total_length = debug_total_length + sub_length;
+        split_sub_list(hap, sub_list, sub_length, hap->snp, overlap_list, R_INF);
+    }
+
+    ///debug_snp_matrix(hap);
+
+    qsort(hap->snp_stat, hap->available_snp, sizeof(SnpStats), cmp_snp_stats);
+
+    ///debug_snp_matrix(hap);
+
+    if(generate_haplotypes(hap))
+    {
+        //print_Haplotype(hap, overlap_list, R_INF);
+        ///if we can phase, we need to exclude the reads of different haplotype
+        int8_t* vector = Get_Result_SNP_Vector((*hap));
+        for (j = 0; j < Get_SNP_Vector_Length((*hap)); j++)
+        {
+            if(vector[j] == 1)
+            {
+                overlap_list->list[j].align_length = 0;
+            }
+        }
+    }
+
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
 
 void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF, 
                         UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read, Graph* g,
                         long long* matched_overlap_0, long long* matched_overlap_1, 
-                        long long* potiental_matched_overlap_0, long long* potiental_matched_overlap_1)
+                        long long* potiental_matched_overlap_0, long long* potiental_matched_overlap_1,
+                        Cigar_record* current_cigar, haplotype_evdience_alloc* hap)
 {
     reverse_complement(g_read->seq, g_read->length);
 
@@ -2581,6 +3882,18 @@ void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF,
             fprintf(stderr, "error length\n");
         }
         
+        /**
+        if(memcmp("m54334_180926_225337/39780640/ccs", Get_NAME((*R_INF), overlap_list->list[0].x_id), 
+    Get_NAME_LENGTH((*R_INF), overlap_list->list[0].x_id)) == 0)
+        {
+            fprintf(stderr, "window_start: %d, window_end: %d, dumy->length: %d, dumy->lengthNT: %d\n", 
+            window_start, window_end, dumy->length, dumy->lengthNT);
+        }
+        **/
+        
+        
+
+        
         ///verify_get_interval(window_start, window_end, overlap_list, dumy);
 
         verify_window(window_start, window_end, overlap_list, dumy, R_INF, g_read->seq);
@@ -2593,15 +3906,58 @@ void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF,
         }
     }
 
-    
     debug_stats(overlap_list, R_INF, g_read, dumy, overlap_read, potiental_matched_overlap_0, potiental_matched_overlap_1);
     
     recalcate_window(overlap_list, R_INF, g_read, dumy, overlap_read);
     
     debug_stats(overlap_list, R_INF, g_read, dumy, overlap_read, matched_overlap_0, matched_overlap_1);
 
+    /**
+        if(memcmp("m54334_180926_225337/39780640/ccs", Get_NAME((*R_INF), overlap_list->list[0].x_id), 
+    Get_NAME_LENGTH((*R_INF), overlap_list->list[0].x_id)) == 0)
+    {
+        fprintf(stderr, "start_pos: %d, end_pos: %d, x_length: %d\n", 
+        g_read->length - 1 - 84, g_read->length - 1, g_read->length);
+
+
+        for (i = 0; i < overlap_list->length; i++)
+        {
+            int j;
+            for(j = 0; j < overlap_list->list[i].w_list_length; j++)
+            {
+                if(overlap_list->list[i].w_list[j].x_start == 13500)
+                {
+                    if(overlap_list->list[i].w_list[j].y_end != -1)
+                    {
+                        fprintf(stderr, "y_name: %.*s, x_start: %d, x_end: %d, y_start: %d, y_end: %d\n", 
+                        Get_NAME_LENGTH((*R_INF), overlap_list->list[i].y_id), Get_NAME((*R_INF), overlap_list->list[i].y_id),
+                        overlap_list->list[i].w_list[j].x_start, overlap_list->list[i].w_list[j].x_end,
+                        overlap_list->list[i].w_list[j].y_start,
+                        overlap_list->list[i].w_list[j].y_end);
+
+                    
+                        int ijk;
+                        for (ijk = 0; ijk < overlap_list->list[i].w_list[j].cigar.length; ijk++)
+                        {
+                            fprintf(stderr, "operation: %d, length: %d\n",
+                            overlap_list->list[i].w_list[j].cigar.C_C[ijk],
+                            overlap_list->list[i].w_list[j].cigar.C_L[ijk]
+                            );
+                        }
+                        
+                    }
+                }
+            }
+        }
+    }
+    **/
+    
+    
+    partition_overlaps(overlap_list, R_INF, g_read, dumy, hap);
+
+    generate_consensus(overlap_list, R_INF, g_read, dumy, g, current_cigar);
+
 
-    generate_consensus(overlap_list, R_INF, g_read, dumy, g);
 
 
 
@@ -2676,6 +4032,41 @@ void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF,
 }
 
 
+void init_Cigar_record(Cigar_record* dummy)
+{
+    dummy->length = 0;
+    dummy->size = 100;
+    dummy->record = (uint32_t*)malloc(sizeof(uint32_t)*dummy->size);
+
+
+    dummy->lost_base_length = 0;
+    dummy->lost_base_size = 100;
+    dummy->lost_base = (char*)malloc(sizeof(char)*dummy->lost_base_size);
+
+
+
+    dummy->current_operation_length = 0;
+    dummy->current_operation = 127;
+}
+
+
+void destory_Cigar_record(Cigar_record* dummy)
+{
+    free(dummy->record);
+    free(dummy->lost_base);
+}
+
+void clear_Cigar_record(Cigar_record* dummy)
+{
+    dummy->new_read_length = 0;
+    dummy->length = 0;
+    dummy->lost_base_length = 0;
+
+    dummy->current_operation_length = 0;
+    dummy->current_operation = 127;
+}
+
+
 void init_Correct_dumy(Correct_dumy* list)
 {
     list->size = 0;
diff --git a/Correct.h b/Correct.h
index e147e32..3ad3f0e 100644
--- a/Correct.h
+++ b/Correct.h
@@ -4,6 +4,7 @@
 #include "Hash_Table.h"
 #include "Levenshtein_distance.h"
 #include "POA.h"
+#include "Process_Read.h"
 
 #define CORRECT_THRESHOLD 0.7
 #define MIN_COVERAGE_THRESHOLD 4
@@ -12,11 +13,324 @@
 #define INSERTION 2
 #define DELETION 3
 
+#define FLAG_THRE 1
+
 #define MAX(x, y) ((x >= y)?x:y)
 #define MIN(x, y) ((x <= y)?x:y)
 #define OVERLAP(x_start, x_end, y_start, y_end) (MIN(x_end, y_end) - MAX(x_start, y_start) + 1) 
 ///#define OVERLAP(x_start, x_end, y_start, y_end) MIN(x_end, y_end) - MAX(x_start, y_start) + 1 
 
+#define Get_MisMatch_Base(RECORD) (s_H[(RECORD>>3)])
+#define Get_Match_Base(RECORD) (s_H[(RECORD&7)])
+
+
+
+typedef struct
+{
+    /**[0-1] bits are type:**/
+    /**[2-31] bits are length**/
+    char current_operation;
+    int current_operation_length;
+    uint32_t* record;
+    uint64_t size;
+    uint64_t length;
+    uint32_t new_read_length;
+
+
+    char* lost_base;
+    uint64_t lost_base_size;
+    uint64_t lost_base_length;
+    
+
+}Cigar_record;
+
+
+typedef struct
+{
+    ////the position of snp in read itself
+    uint32_t site;
+    ////the overlapID
+    uint32_t overlapID;
+    ////the position of snp in that overlap
+    uint32_t overlapSite;
+    ///there are several types: 0: equal to read 1: not equal to read, but it is a mismatch 2: is a gap
+    uint8_t type;
+    ///misbase
+    char misBase;
+}haplotype_evdience;
+
+
+
+typedef struct
+{
+    ///the id of this snp
+    uint32_t id;
+    uint32_t overlap_num;
+    uint32_t occ_0;
+    uint32_t occ_1;
+    uint32_t occ_2;
+    int score;
+    ////the position of snp in read itself
+    uint32_t site;
+}
+SnpStats;
+
+
+#define Get_SNP_Martix_Size(matrix) (matrix.snp * matrix.overlap)
+#define Get_SNP_Vector(matrix, i) (matrix.snp_matrix + matrix.overlap * i)
+#define Get_SNP_Vector_Length(matrix) (matrix.overlap)
+#define Get_Result_SNP_Vector(matrix) (matrix.snp_matrix + matrix.overlap*matrix.snp)
+
+typedef struct
+{
+    haplotype_evdience* list;
+    uint32_t sub_list_start;
+    uint32_t sub_list_length;
+    uint32_t length;
+    uint32_t size;
+    
+    uint32_t flag[WINDOW];
+
+    uint32_t available_snp;
+    uint32_t core_snp;
+    uint32_t snp;
+    uint32_t overlap;
+    int8_t* snp_matrix;
+    uint32_t snp_matrix_size;
+    SnpStats* snp_stat; 
+    SnpStats result_stat;
+    uint32_t snp_stat_size;
+}
+haplotype_evdience_alloc;
+
+inline int filter_snp(int x, int y, int total)
+{
+    double available;
+
+    if(x <= y)
+    {
+        available = x;
+    }
+    else
+    {
+        available = y;
+    }
+    double threshold = 0.30;
+    available = available/((double)(total));
+    if(available <= threshold && available < 6)
+    {
+        return 0;
+    }
+    return 1;
+}
+
+
+inline int filter_one_snp(int occ_0, int occ_1, int total)
+{
+    
+
+    double available;
+
+    if(occ_0 <= occ_1)
+    {
+        available = occ_0;
+    }
+    else
+    {
+        available = occ_1;
+    }
+    double threshold = 0.35;
+    available = available/((double)(total));
+
+    
+
+    if(available < threshold || occ_0 < MIN_COVERAGE_THRESHOLD + 1 || total < 10)
+    {
+        return 0;
+    }
+    return 1;
+}
+
+inline void InsertSNPVector(haplotype_evdience_alloc* h, haplotype_evdience* sub_list, long long sub_length, char misBase)
+{
+    if(sub_length <= 0)
+        return;
+    long long i = 0;
+    h->snp_stat[h->available_snp].id = h->available_snp;
+    h->snp_stat[h->available_snp].occ_0 = 0;
+    h->snp_stat[h->available_snp].occ_1 = 0;
+    h->snp_stat[h->available_snp].occ_2 = 0;
+    h->snp_stat[h->available_snp].overlap_num = 0;
+    
+    h->snp_stat[h->available_snp].site = sub_list[0].site;
+
+    int8_t* vector = Get_SNP_Vector((*h), h->available_snp);
+    for (i = 0; i < sub_length; i++)
+    {
+        if(sub_list[i].type == 0)
+        {
+            vector[sub_list[i].overlapID] = 0;
+            h->snp_stat[h->available_snp].occ_0++;
+        }
+        else if(sub_list[i].type == 1 && sub_list[i].misBase == misBase)
+        {
+            vector[sub_list[i].overlapID] = 1;
+            h->snp_stat[h->available_snp].occ_1++;
+        }
+        else
+        {
+            vector[sub_list[i].overlapID] = 2;
+            h->snp_stat[h->available_snp].occ_2++;
+        }
+        h->snp_stat[h->available_snp].overlap_num++;
+    }
+
+
+    int new_occ_0 = h->snp_stat[h->available_snp].occ_0 + 1;
+    int new_occ_1 = h->snp_stat[h->available_snp].occ_1;
+
+    if(filter_snp(new_occ_0, new_occ_1, new_occ_0 + new_occ_1) == 0)
+    {
+        h->snp_stat[h->available_snp].score = -1;
+    }
+    else
+    {
+        h->core_snp++;
+
+
+        double consensus = new_occ_0 + new_occ_1 - abs(new_occ_0 - new_occ_1);
+
+        consensus = consensus /((double)(new_occ_0 + new_occ_1));
+
+        ///50% vs 50%
+        if(new_occ_0 == new_occ_1)
+        {
+            consensus = consensus + 0.25;
+        }
+        else if(consensus >= 0.8)
+        {
+            consensus = consensus + 0.2;
+        }
+        else if(consensus >= 0.6)
+        {
+            consensus = consensus + 0.15;
+        }
+        else if(consensus >= 0.4)
+        {
+            consensus = consensus + 0.1;
+        }
+        else if(consensus >= 0.2)
+        {
+            consensus = consensus + 0.05;
+        }
+        
+
+        consensus= consensus*((double)(new_occ_0 + new_occ_1));
+
+        h->snp_stat[h->available_snp].score = consensus;
+    }
+
+
+    
+    h->available_snp++;
+}
+
+inline void SetSnpMatrix(haplotype_evdience_alloc* h, long long snp_num, long long overlap_num)
+{
+    long long new_size = (snp_num + 1)* overlap_num;
+
+    if(h->snp_matrix_size < new_size)
+    {
+        h->snp_matrix_size = new_size;
+        h->snp_matrix = (int8_t*)realloc(h->snp_matrix, h->snp_matrix_size);
+    }
+
+    if(h->snp_stat_size < snp_num)
+    {
+        h->snp_stat_size = snp_num;
+        h->snp_stat = (SnpStats*)realloc(h->snp_stat, h->snp_stat_size * sizeof(SnpStats));
+    }
+
+    ///h->snp may be different with the number of snp vector
+    ///since some snps have been filtered
+    h->snp = snp_num;
+    h->overlap = overlap_num;
+    h->available_snp = 0;
+    h->core_snp = 0;
+
+    memset(h->snp_matrix, -1, h->snp * h->overlap);
+
+    
+}
+
+inline void InitHaplotypeEvdience(haplotype_evdience_alloc* h)
+{
+    h->snp = 0;
+    h->available_snp = 0;
+    h->overlap = 0;
+    h->snp_matrix_size = 0;
+    h->snp_stat_size = 0;
+    h->snp_matrix = NULL;
+    h->snp_stat = NULL;
+
+
+    h->sub_list_start = 0;
+    h->sub_list_length = 0;
+    h->length = 0;
+    h->size = 100;
+    h->list = (haplotype_evdience*)calloc(h->size, sizeof(haplotype_evdience));
+    memset(h->flag, 0, WINDOW * sizeof(uint32_t));
+}
+
+inline void StarSubListHaplotypeEvdience(haplotype_evdience_alloc* h)
+{
+    h->sub_list_start = h->length;
+}
+
+inline void EndSubListHaplotypeEvdience(haplotype_evdience_alloc* h)
+{
+    h->sub_list_length = h->length - h->sub_list_start;
+}
+
+
+inline void destoryHaplotypeEvdience(haplotype_evdience_alloc* h)
+{
+    free(h->list);
+    free(h->snp_stat);
+    free(h->snp_matrix);
+}
+
+inline void ResizeInitHaplotypeEvdience(haplotype_evdience_alloc* h)
+{
+    h->snp = 0;
+    h->length = 0;
+    h->sub_list_start = 0;
+    h->sub_list_length = 0;
+    memset(h->flag, 0, WINDOW * sizeof(uint32_t));
+}
+
+inline void RsetInitHaplotypeEvdienceFlag(haplotype_evdience_alloc* h)
+{
+    memset(h->flag, 0, WINDOW * sizeof(uint32_t));
+}
+
+inline void addHaplotypeEvdience(haplotype_evdience_alloc* h, haplotype_evdience* ev)
+{
+    uint32_t new_length = h->length + 1;
+    if(new_length > h->size)
+    {
+        h->size = h->size * 2;
+        if(h->size < new_length)
+        {
+            h->size = new_length;
+        }
+        h->list = (haplotype_evdience*)realloc(h->list, sizeof(haplotype_evdience)*h->size);
+    }
+
+    h->list[h->length] = (*ev);
+    h->length++;
+}
+
 typedef struct
 {
     char* corrected_read;
@@ -41,7 +355,8 @@ typedef struct
 void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF, 
                         UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read, Graph* g,
                         long long* matched_overlap_0, long long* matched_overlap_1, 
-                        long long* potiental_matched_overlap_0, long long* potiental_matched_overlap_1);
+                        long long* potiental_matched_overlap_0, long long* potiental_matched_overlap_1,
+                        Cigar_record* current_cigar, haplotype_evdience_alloc* hap);
 void init_Correct_dumy(Correct_dumy* list);
 void destory_Correct_dumy(Correct_dumy* list);
 void clear_Correct_dumy(Correct_dumy* list, overlap_region_alloc* overlap_list);
@@ -51,6 +366,192 @@ void pre_filter_by_nearby_single(k_mer_pos* new_n_list, k_mer_pos* old_n_list, u
 All_reads* R_INF, Correct_dumy* dumy, uint64_t* new_n_length);
 void get_seq_from_Graph(Graph* backbone, Correct_dumy* dumy);
 
+void init_Cigar_record(Cigar_record* dummy);
+void destory_Cigar_record(Cigar_record* dummy);
+void clear_Cigar_record(Cigar_record* dummy);
+
+
+
+inline void add_new_cell_to_cigar_record(Cigar_record* dummy, uint32_t len, uint32_t type)
+{
+    uint32_t tmp;
+    tmp = len;
+    tmp = tmp << 2;
+    tmp = tmp | type;
+
+    dummy->length++;
+
+
+    if(dummy->length > dummy->size)
+    {
+        dummy->size = dummy->size * 2;
+        dummy->record = (uint32_t*)realloc(dummy->record, dummy->size*sizeof(uint32_t));
+    }
+
+    dummy->record[dummy->length - 1] = tmp;
+}
+
+inline void add_existing_cell_to_cigar_record(Cigar_record* dummy, uint32_t len, uint32_t type)
+{
+    uint32_t tmp;
+
+    tmp = dummy->record[dummy->length - 1] >> 2;
+    tmp = tmp + len;
+    tmp = tmp << 2;
+    tmp = tmp | type;
+    dummy->record[dummy->length - 1] = tmp;
+}
+
+
+inline void add_new_cell_to_cigar_record_with_different_base(Cigar_record* dummy, uint32_t len, uint32_t type, char* seq)
+{
+    uint32_t tmp;
+    tmp = len;
+    tmp = tmp << 2;
+    tmp = tmp | type;
+   
+
+    dummy->length++;
+
+    if(dummy->length > dummy->size)
+    {
+        dummy->size = dummy->size * 2;
+        dummy->record = (uint32_t*)realloc(dummy->record, dummy->size*sizeof(uint32_t));
+    }
+
+    dummy->record[dummy->length - 1] = tmp;
+
+    
+
+    if (dummy->lost_base_length + len> dummy->lost_base_size)
+    {
+        dummy->lost_base_size = (dummy->lost_base_length + len) * 2;
+        dummy->lost_base = (char*)realloc(dummy->lost_base, dummy->lost_base_size*sizeof(char));
+    }
+
+    int i = 0;
+    for (i = 0; i < len; i++, dummy->lost_base_length++)
+    {
+        dummy->lost_base[dummy->lost_base_length] = seq[i];
+    }
+}
+
+inline void add_existing_cell_to_cigar_record_with_different_base(Cigar_record* dummy, uint32_t len, uint32_t type, char* seq)
+{
+    uint32_t tmp;
+
+    tmp = dummy->record[dummy->length - 1] >> 2;
+    tmp = tmp + len;
+    tmp = tmp << 2;
+    tmp = tmp | type;
+    dummy->record[dummy->length - 1] = tmp;
+
+    if (dummy->lost_base_length + len> dummy->lost_base_size)
+    {
+        dummy->lost_base_size = (dummy->lost_base_length + len) * 2;
+        dummy->lost_base = (char*)realloc(dummy->lost_base, dummy->lost_base_size*sizeof(char));
+    }
+
+    int i = 0;
+    for (i = 0; i < len; i++, dummy->lost_base_length++)
+    {
+        dummy->lost_base[dummy->lost_base_length] = seq[i];
+    }
+}
+
+
+/***
+ type:
+ 0. match
+ 1. mismatch
+ 2. insertion
+ 3. deletion
+ ***/
+inline void add_cigar_record(char* seq, uint32_t len, Cigar_record* dummy, uint32_t type)
+{
+
+    
+    uint32_t tmp;
+    
+    if(type == 0)///match
+    {
+        ///add to existing cell, just increase length
+        if(dummy->current_operation == type)
+        {
+            add_existing_cell_to_cigar_record(dummy, len, type);
+        }
+        else  ///add to new cell
+        {
+            add_new_cell_to_cigar_record(dummy, len, type);
+        }
+
+        dummy->new_read_length += len;
+    }
+    else if(type == 1)///mismatch
+    {
+        ///add to existing cell, just increase length
+        ///and add different bases
+        if(dummy->current_operation == type)
+        {
+            add_existing_cell_to_cigar_record_with_different_base(dummy, len, type, seq);
+        }
+        else
+        {
+            add_new_cell_to_cigar_record_with_different_base(dummy, len, type, seq);
+        }
+
+        dummy->new_read_length += len;
+    }
+    else if(type == 3)///deletion, the bases in previous read will be removed
+    {
+        ///add to existing cell, just increase length
+        ///and add different bases
+        if(dummy->current_operation == type)
+        {
+            add_existing_cell_to_cigar_record_with_different_base(dummy, len, type, seq);
+        }
+        else
+        {
+            add_new_cell_to_cigar_record_with_different_base(dummy, len, type, seq);
+        }
+    }
+    else if(type == 2)///insertion
+    {
+        /**
+        ///add to existing cell, just increase length
+        if(dummy->current_operation == type)
+        {
+            add_existing_cell_to_cigar_record(dummy, len, type);
+        }
+        else  ///add to new cell
+        {
+            add_new_cell_to_cigar_record(dummy, len, type);
+        }
+        **/
+       ///add to existing cell, just increase length
+        ///and add different bases
+        if(dummy->current_operation == type)
+        {
+            add_existing_cell_to_cigar_record_with_different_base(dummy, len, type, seq);
+        }
+        else
+        {
+            add_new_cell_to_cigar_record_with_different_base(dummy, len, type, seq);
+        }
+
+        dummy->new_read_length += len;
+    }
+    
+    dummy->current_operation = type;
+
+
+
+
+
+    
+}
+
+
 /**********************for prefilter************************ */
 void destory_k_mer_pos_list_alloc_prefilter(k_mer_pos_list_alloc* list);
 void append_k_mer_pos_list_alloc_prefilter(k_mer_pos_list_alloc* list, k_mer_pos* n_list, uint64_t n_length, 
diff --git a/Hash_Table.cpp b/Hash_Table.cpp
index 61eef07..9209e8c 100644
--- a/Hash_Table.cpp
+++ b/Hash_Table.cpp
@@ -658,7 +658,8 @@ uint64_t readID, uint64_t readLength, All_reads* R_INF)
 
 
 
-void append_window_list(overlap_region* region, uint64_t x_start, uint64_t x_end, int y_start, int y_end, int error)
+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)
 {
     
     long long length = region->x_pos_e - region->x_pos_s + 1;
@@ -678,6 +679,8 @@ void append_window_list(overlap_region* region, uint64_t x_start, uint64_t x_end
     region->w_list[region->w_list_length].y_end = y_end;
     region->w_list[region->w_list_length].error = error;
     region->w_list[region->w_list_length].cigar.length = -1;
+    region->w_list[region->w_list_length].extra_begin = extra_begin;
+    region->w_list[region->w_list_length].extra_end = extra_end;
     region->w_list_length++;
 }
 
diff --git a/Hash_Table.h b/Hash_Table.h
index e1fbb61..8520b97 100644
--- a/Hash_Table.h
+++ b/Hash_Table.h
@@ -82,10 +82,13 @@ typedef struct
 
 typedef struct
 {
+  ///the begining and end of a window, instead of the whole overlap
   uint64_t x_start;
   uint64_t x_end;
   int y_end;
   int y_start;
+  int extra_begin;
+  int extra_end;
   ///int y_pre_start;
   ///error小于等于0都要重新算
   int error;
@@ -95,6 +98,7 @@ typedef struct
 typedef struct
 {
     uint64_t x_id;
+    ///the begining and end of the whole overlap
     uint64_t x_pos_s;
     uint64_t x_pos_e;
     uint64_t x_pos_strand;
@@ -450,7 +454,8 @@ void destory_overlap_region_alloc(overlap_region_alloc* list);
 void append_overlap_region_alloc(overlap_region_alloc* list, overlap_region* tmp, All_reads* R_INF);
 void calculate_overlap_region(Candidates_list* candidates, overlap_region_alloc* overlap_list, 
 uint64_t readID, uint64_t readLength, All_reads* R_INF);
-void append_window_list(overlap_region* region, uint64_t x_start, uint64_t x_end, int y_start, int y_end, int error);
+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);
 
 
 
diff --git a/POA.cpp b/POA.cpp
index b7ed636..508a900 100644
--- a/POA.cpp
+++ b/POA.cpp
@@ -256,8 +256,6 @@ void addmatchedSeqToGraph(Graph* backbone, long long currentNodeID, char* x_stri
     int last_operation = -1;
 
     
-    
-    
     ///note that node 0 is the start node
     ///0 is match, 1 is mismatch, 2 is up, 3 is left
     ///2是x缺字符（y多字符），而3是y缺字符（x多字符）
diff --git a/Process_Read.cpp b/Process_Read.cpp
index c08f9e6..957fec0 100644
--- a/Process_Read.cpp
+++ b/Process_Read.cpp
@@ -186,6 +186,9 @@ int load_All_reads(All_reads* r, char* read_file_name)
 	r->read_length = (uint64_t*)malloc(sizeof(uint64_t)*r->total_reads);
 	fread(r->read_length, sizeof(uint64_t), r->total_reads, fp);
 
+	r->read_size = (uint64_t*)malloc(sizeof(uint64_t)*r->total_reads);
+	memcpy (r->read_size, r->read_length, sizeof(uint64_t)*r->total_reads);
+
 	/**********should remove**********/
 	///r->read = (uint8_t*)malloc(sizeof(uint8_t)*(r->total_reads_bases/4 + r->total_reads + 5));
 	///fread(r->read, sizeof(uint8_t), (r->total_reads_bases/4 + r->total_reads + 5), fp);
@@ -205,6 +208,19 @@ int load_All_reads(All_reads* r, char* read_file_name)
 	fread(r->name_index, sizeof(uint64_t), r->name_index_size, fp);
 
 
+	r->cigars = (Compressed_Cigar_record*)malloc(sizeof(Compressed_Cigar_record)*r->total_reads);
+	for (i = 0; i < r->total_reads; i++)
+	{
+		r->cigars[i].size = 0;
+		r->cigars[i].length = 0;
+		r->cigars[i].record = NULL;
+
+		r->cigars[i].lost_base_size = 0;
+		r->cigars[i].lost_base_length = 0;
+		r->cigars[i].lost_base = NULL;
+	}
+
+
     free(index_name);    
     fclose(fp);
     fprintf(stdout, "Reads has been loaded.\n");
@@ -246,6 +262,9 @@ inline void insert_read(All_reads* r, kstring_t* read, kstring_t* name)
 void malloc_All_reads(All_reads* r)
 {
 
+	r->read_size = (uint64_t*)malloc(sizeof(uint64_t)*r->total_reads);
+	memcpy (r->read_size, r->read_length, sizeof(uint64_t)*r->total_reads);
+
 	///必须加r->total_reads
 	/**********should remove**********/
 	///r->read = (uint8_t*)malloc(sizeof(uint8_t)*(r->total_reads_bases/4 + r->total_reads + 5));
@@ -257,7 +276,17 @@ void malloc_All_reads(All_reads* r)
 		r->read_sperate[i] = (uint8_t*)malloc(sizeof(uint8_t)*(r->read_length[i]/4+1));
 	}
 
+	r->cigars = (Compressed_Cigar_record*)malloc(sizeof(Compressed_Cigar_record)*r->total_reads);
+	for (i = 0; i < r->total_reads; i++)
+	{
+		r->cigars[i].size = 0;
+		r->cigars[i].length = 0;
+		r->cigars[i].record = NULL;
 
+		r->cigars[i].lost_base_size = 0;
+		r->cigars[i].lost_base_length = 0;
+		r->cigars[i].lost_base = NULL;
+	}
 
 
 	r->name = (char*)malloc(sizeof(char)*r->total_name_length);
@@ -298,6 +327,119 @@ void init_UC_Read(UC_Read* r)
 }
 
 
+void recover_UC_Read_sub_region_begin_end
+(char* r, long long start_pos, long long length, uint8_t strand, All_reads* R_INF, long long ID, int extra_begin, int extra_end)
+{
+
+
+
+	long long readLen = Get_READ_LENGTH((*R_INF), ID);
+	uint8_t* src = Get_READ((*R_INF), ID);
+
+	long long i;
+	long long copyLen;
+	long long end_pos = start_pos + length - 1;
+
+
+	
+
+	if (strand == 0)
+	{
+		
+		i = start_pos;
+		copyLen = 0;
+
+
+
+		long long initLen = start_pos % 4;
+
+		if (initLen != 0)
+		{
+			memcpy(r, bit_t_seq_table[src[i>>2]] + initLen, 4 - initLen);
+			copyLen = copyLen + 4 - initLen;
+			i = i + copyLen;
+		}
+		
+		
+		while (copyLen < length)
+		{
+			memcpy(r+copyLen, bit_t_seq_table[src[i>>2]], 4);
+			copyLen = copyLen + 4;
+			i = i + 4;
+		}
+
+		
+		if (R_INF->N_site[ID])
+		{
+			for (i = 1; i <= R_INF->N_site[ID][0]; i++)
+			{
+				if (R_INF->N_site[ID][i] >= start_pos && R_INF->N_site[ID][i] <= end_pos)
+				{
+					r[R_INF->N_site[ID][i] - start_pos] = 'N';
+				}
+				else if(R_INF->N_site[ID][i] > end_pos)
+				{
+					break;
+				}
+			}
+		}
+		
+		
+	}
+	else
+	{
+		
+		start_pos = readLen - start_pos - 1;
+		end_pos = readLen - end_pos - 1;
+
+
+
+		///start_pos > end_pos
+		i = start_pos;
+		copyLen = 0;
+		long long initLen = (start_pos + 1) % 4;
+
+		if (initLen != 0)
+		{
+			memcpy(r, bit_t_seq_table_rc[src[i>>2]] + 4 - initLen, initLen);
+			copyLen = copyLen + initLen;
+			i = i - initLen;
+		}
+
+		while (copyLen < length)
+		{
+			memcpy(r+copyLen, bit_t_seq_table_rc[src[i>>2]], 4);
+			copyLen = copyLen + 4;
+			i = i - 4;
+		}
+
+		if (R_INF->N_site[ID])
+		{
+			long long offset = readLen - start_pos - 1;
+
+			for (i = 1; i <= R_INF->N_site[ID][0]; i++)
+			{
+
+				if (R_INF->N_site[ID][i] >= end_pos && R_INF->N_site[ID][i] <= start_pos)
+				{
+					r[readLen - R_INF->N_site[ID][i] - 1 - offset] = 'N';
+				}
+				else if(R_INF->N_site[ID][i] > start_pos)
+				{
+					break;
+				}
+			}
+		}
+		
+
+		
+	}
+		
+}
+
+
+
+
 void recover_UC_Read_sub_region(char* r, long long start_pos, long long length, uint8_t strand, All_reads* R_INF, long long ID)
 {
 
diff --git a/Process_Read.h b/Process_Read.h
index 4981e04..b303f33 100644
--- a/Process_Read.h
+++ b/Process_Read.h
@@ -48,8 +48,8 @@ static uint8_t seq_nt6_table[256] = {
 
 static char bit_t_seq_table[256][4] = {0};
 static char bit_t_seq_table_rc[256][4] = {0};
-static char s_H[4] = {'A', 'C', 'G', 'T'};
-static char rc_Table[4] = {'T', 'G', 'C', 'A'};
+static char s_H[5] = {'A', 'C', 'G', 'T', 'N'};
+static char rc_Table[5] = {'T', 'G', 'C', 'A', 'N'};
 
 #define RC_CHAR(x) rc_Table[seq_nt6_table[(uint8_t)x]]
 
@@ -58,6 +58,19 @@ void init_kseq(char* file);
 void destory_kseq();
 int get_read(kseq_t *s);
 
+typedef struct
+{
+    /**[0-1] bits are type:**/
+    /**[2-31] bits are length**/
+    uint32_t* record;
+    uint32_t length;
+	uint32_t size;
+
+    char* lost_base;
+    uint32_t lost_base_length;
+	uint32_t lost_base_size;
+	uint32_t new_length;
+}Compressed_Cigar_record;
 
 
 typedef struct
@@ -69,6 +82,7 @@ typedef struct
 
 	uint8_t** read_sperate;
 	uint64_t* read_length;
+	uint64_t* read_size;
 
 	///seq start pos in uint8_t* read
 	///do not need it
@@ -83,6 +97,8 @@ typedef struct
 	uint64_t total_reads_bases;
 	uint64_t total_name_length;
 
+	Compressed_Cigar_record* cigars; 
+
 } All_reads;
 
 extern All_reads R_INF;
diff --git a/debug_assembly.sh b/debug_assembly.sh
new file mode 100644
index 0000000..506f1a8
--- /dev/null
+++ b/debug_assembly.sh
@@ -0,0 +1,17 @@
+#!/bin/bash
+# My first script
+if [ $# -eq 1 ]
+then
+echo "../minimap2/minimap2 -ax asm20 -t 32 ../minimap2/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz "$1" >"$1".sam"
+../minimap2/minimap2 -ax asm20 -t 32 ../minimap2/Homo_sapiens.GRCh38.dna.primary_assembly.fa $1 >$1.sam
+echo "samtools view -Sb "$1".sam >"$1".bam"
+samtools view -Sb $1.sam >$1.bam
+echo "samtools sort "$1".bam sort_"$1
+samtools sort $1.bam sort_$1
+echo "rm sort_"$1".bam.bai"
+rm sort_$1.bam.bai
+echo "samtools index sort_"$1".bam"
+samtools index sort_$1.bam
+else
+echo "debug_assembly.sh intput.fa"
+fi
diff --git a/main.cpp b/main.cpp
index 11ea66f..97333d5 100644
--- a/main.cpp
+++ b/main.cpp
@@ -206,6 +206,7 @@ int main(int argc, char *argv[])
 
     return 1;
     **/
+   fprintf(stdout, "Will perform %d round of error correction...\n", number_of_round);
     
    fprintf(stdout, "defined k_mer_min_freq by user: %d\n", k_mer_min_freq);
    fprintf(stdout, "defined k_mer_max_freq by user: %d\n", k_mer_max_freq);
@@ -213,25 +214,7 @@ int main(int argc, char *argv[])
 
     fprintf(stdout, "k-mer length: %d\n",k_mer_length);
 
-    if (load_index_from_disk && load_pre_cauculated_index())
-    {
-        ;
-    }
-    else
-    {
-        Counting_multiple_thr();
-
-        Build_hash_table_multiple_thr();
-    }
-
-    fprintf(stdout, "k_mer_min_freq in hashtable: %d\n", k_mer_min_freq);
-    fprintf(stdout, "k_mer_max_freq in hashtable: %d\n", k_mer_max_freq);
-
-    
-    
-    ///verify_Position_hash_table();
-
-    Overlap_calculate_multipe_thr();
+    Correct_Reads(number_of_round);
 
 
 
diff --git a/sort_fasta.sh b/sort_fasta.sh
new file mode 100644
index 0000000..15d9e6e
--- /dev/null
+++ b/sort_fasta.sh
@@ -0,0 +1,10 @@
+#!/bin/bash
+# My first script
+if [ $# -eq 1 ]
+then
+echo "input file is: "$1
+echo $1" | paste - - | sort -k1,1 -t \" \" | tr \"\t\" \"\n\" > sort_"$1
+cat $1 | paste - - | sort -k1,1 -t " " | tr "\t" "\n" > sort_$1
+else
+echo "./sort_fasta.sh file_name"
+fi