Finish multiple thread counting

Assembly.cpp

@@ -4,11 +4,15 @@
 #include <zlib.h>
 #include "Process_Read.h"
 #include "CommandLines.h"
+#include "kmer.h"
+#include "Hash_Table.h"
+
+Total_Count_Table TCB;
 
 
 
-
-void Counting()
+/********************************for debug***************************************/
+void Verify_Counting()
 {
 
 
@@ -17,46 +21,87 @@ void Counting()
 
     long long read_number = 0;
 
+    HPC_seq HPC_read;
+    uint64_t code;
+    Hash_code k_code;
+    int avalible_k = 0;
+
+    fprintf(stdout, "Start Verifying ...\n");
+
     while (get_read(seq))
     {
 
-        fprintf(stderr,"@%s\n",seq->name.s);
-        fprintf(stderr,"%s\n",seq->seq.s);
-        fprintf(stderr,"+\n");
-        fprintf(stderr,"%s\n",seq->qual.s);
+        init_HPC_seq(&HPC_read, seq->seq.s, seq->seq.l);
+        init_Hash_code(&k_code);
+
+        avalible_k = 0;
+
+        while ((code = get_HPC_code(&HPC_read)) != 6)
+        {
+            if(code < 4)
+            {
+                k_mer_append(&k_code,code,k_mer_length);
+                avalible_k++;
+                if (avalible_k>=k_mer_length)
+                {
+                    if(verify_Total_Count_Table(&TCB, &k_code, k_mer_length) == -1)
+                    {
+                        fprintf(stderr, "ERROR when subtracting!\n");
+                    }
+                }
+                
+            }
+            else
+            {
+                avalible_k = 0;
+                init_Hash_code(&k_code);
+            }
+            
+        }
+
+
+
 
         read_number++;
     }
     
     fprintf(stdout, "read_number: %lld\n",read_number);
 
+    fprintf(stdout, "Start Traversing ...\n");
+
+    Traverse_Total_Count_Table(&TCB);
+
+    fprintf(stdout, "Finish Traversing!\n");
 
 
 }
-
-
-
-
-
-
-
-
-
-
-
+/********************************for debug***************************************/
 
 
 
 
 void* Perform_Counting(void* arg)
 {
+    int thr_ID = *((int*)arg);
+    int i = 0;
+    HPC_seq HPC_read;
+
     R_buffer_block curr_sub_block;
 
     init_R_buffer_block(&curr_sub_block);
 
     long long read_number = 0;
-    
+    long long k_mer_number = 0 ;
+
 	int file_flag = 1;
+
+    uint64_t code;
+
+    Hash_code k_code;
+
+    int avalible_k = 0;
+
+
 	while (file_flag != 0)
 	{
 
@@ -64,10 +109,46 @@ void* Perform_Counting(void* arg)
 
         read_number = read_number + curr_sub_block.num; 
 
+        for (i = 0; i < curr_sub_block.num; i++)
+        {
+
+            init_HPC_seq(&HPC_read, curr_sub_block.read[i].seq.s, curr_sub_block.read[i].seq.l);
+            init_Hash_code(&k_code);
+
+            avalible_k = 0;
+
+            while ((code = get_HPC_code(&HPC_read)) != 6)
+            {
+                if(code < 4)
+                {
+                    k_mer_append(&k_code,code,k_mer_length);
+                    avalible_k++;
+                    if (avalible_k>=k_mer_length)
+                    {
+                        ///插入
+                        insert_Total_Count_Table(&TCB, &k_code, k_mer_length);
+                        k_mer_number++;
+                    }
+                    
+                }
+                else
+                {
+                    avalible_k = 0;
+                    init_Hash_code(&k_code);
+                }
+                
+            }
+            
+
+        }
+        
+
     }
 
-    fprintf(stdout, "#########read_number: %lld\n",read_number);
-    fflush(stdout);
+    destory_R_buffer_block(&curr_sub_block);
+    free(arg);
+
+    fprintf(stdout, "thr_ID: %d, read_number: %lld, k_mer_number: %lld\n",thr_ID, read_number, k_mer_number);
 
 }
 
@@ -89,6 +170,15 @@ void* Perform_Counting(void* arg)
 
 void Counting_multiple_thr()
 {
+
+    init_Total_Count_Table(k_mer_length, &TCB);
+
+    fprintf(stdout, "TCB.prefix_bits: %d\n", TCB.prefix_bits);
+    fprintf(stdout, "TCB.suffix_bits: %d\n", TCB.suffix_bits);
+    fprintf(stdout, "TCB.size: %d\n", TCB.size);
+
+
+
     pthread_t inputReadsHandle;
     
     init_R_buffer(thread_num);
@@ -120,6 +210,11 @@ void Counting_multiple_thr()
 
     free(_r_threads);
 
+    destory_R_buffer();
+
+
+    destory_Total_Count_Table(&TCB);
+    
 
 }
 
@@ -134,6 +229,14 @@ void Counting_multiple_thr()
 
 
 
+
+
+
+
+
+
+
+
 
 
 

Assembly.h

@@ -3,7 +3,9 @@
 
 
 
-void Counting();
 void Counting_multiple_thr();
+/********************************for debug***************************************/
+void Verify_Counting();
+/********************************for debug***************************************/
 
 #endif

CommandLines.cpp

@@ -6,6 +6,7 @@
 char* read_file_name = NULL;
 char* output_file_name = NULL;
 int thread_num = 1;
+int k_mer_length = 40;
 
 void Print_H()
 {
@@ -29,11 +30,12 @@ int CommandLine_process (int argc, char *argv[])
     ketopt_t opt = KETOPT_INIT;
 
     int i, c;
-    while ((c = ketopt(&opt, argc, argv, 1, "ht:o:q:", longopts)) >= 0) {
+    while ((c = ketopt(&opt, argc, argv, 1, "ht:o:q:k:", 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 == 'k') k_mer_length = atoi(opt.arg); 
         else if (c == '?') printf("unknown opt: -%c\n", opt.opt? opt.opt : ':');
         else if (c == ':') printf("missing arg: -%c\n", opt.opt? opt.opt : ':');
     }

CommandLines.h

@@ -7,6 +7,8 @@
 extern char* read_file_name;
 extern char* output_file_name;
 extern int thread_num;
+extern int k_mer_length;
+
 
 int CommandLine_process (int argc, char *argv[]);
 

Hash_Table.cpp

@@ -0,0 +1,114 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <stdint.h>
+#include "Hash_Table.h"
+
+
+void init_Count_Table(Count_Table** table)
+{
+    *table = kh_init(COUNT64);
+}
+
+void init_Total_Count_Table(int k, Total_Count_Table* TCB)
+{
+    if(k>64)
+    {
+        fprintf(stdout, "k-mer is too long. The length of k-mer must <= 64.");
+        fflush(stdout);
+        exit(0);
+    }
+
+    int total_bits = k * 2; 
+    TCB->prefix_bits = PREFIX_BITS;
+    TCB->suffix_bits = total_bits - TCB->prefix_bits;
+    if (TCB->suffix_bits > MAX_SUFFIX_BITS)
+    {
+        TCB->suffix_bits = MAX_SUFFIX_BITS;
+        TCB->prefix_bits = total_bits - TCB->suffix_bits;
+    }
+    ///TCB->suffix_mode = (1ULL<<TCB->suffix_bits) - 1;
+    ///这个右移是安全的，因为TCB->suffix_bits不可能是0
+    TCB->suffix_mode = ALL >> (64 - TCB->suffix_bits);
+    
+    ///number of small hash table
+    TCB->size = (1ULL<<TCB->prefix_bits);
+    TCB->sub_h = (Count_Table**)malloc(sizeof(Count_Table*)*TCB->size);
+    TCB->sub_h_lock = (Hash_table_spin_lock*)malloc(sizeof(Hash_table_spin_lock)*TCB->size);
+
+    int i = 0;
+    for (i = 0; i < TCB->size; i++)
+    {
+        init_Count_Table(&(TCB->sub_h[i]));
+        TCB->sub_h_lock[i].lock = 0;
+    }
+
+
+
+}
+
+void destory_Total_Count_Table(Total_Count_Table* TCB)
+{
+    int i;
+    for (i = 0; i < TCB->size; i++)
+    {
+        kh_destroy(COUNT64, TCB->sub_h[i]);
+    }
+    free(TCB->sub_h);
+    free(TCB->sub_h_lock);
+}
+
+
+
+
+
+/********************************for debug***************************************/
+void test_COUNT64()
+{
+
+    uint64_t sb[7] = {5000000000, 6000000000, 7000000000, 8000000000, 9000000000, 10000000000 ,5000000000};
+    Count_Table* h;
+    ///构建哈希表
+    init_Count_Table(&h);
+
+    int i;
+    khint_t k;  ///这就是个迭代器
+    int absent;
+
+    for (i = 0; i < 7; i++)
+    {
+        ///将5作为key插入到哈希表COUNT64中
+        ///absent为0代表哈希表里面已经有这个key了
+        k = kh_put(COUNT64, h, sb[i], &absent);
+        ///代表插入了一个哈希表中没有的元素
+        ///则直接插入
+        if (absent)
+        {
+            kh_value(h, k) = i;
+        }
+        else   ///哈希表中已有的元素
+        {
+            kh_value(h, k) = i;
+        }
+    }
+
+    for (i = 0; i < 7; i++)
+    {
+        ///查询哈希表，key为k
+        k = kh_get(COUNT64, h, sb[i]);
+
+        if (k != kh_end(h))
+        {
+            fprintf(stderr, "value: %d\n", kh_value(h, k));
+        }
+        else
+        {
+            fprintf(stderr, "not found!\n");
+        }
+    }
+    
+
+    kh_destroy(COUNT64, h);
+    
+
+}
+/********************************for debug***************************************/

Hash_Table.h

@@ -0,0 +1,189 @@
+#ifndef __HASHTABLE__
+#define __HASHTABLE__
+#include "khash.h"
+#include "kmer.h"
+
+KHASH_MAP_INIT_INT64(COUNT64, int) 
+
+typedef khash_t(COUNT64) Count_Table;
+
+#define PREFIX_BITS 16
+#define MAX_SUFFIX_BITS 64
+
+typedef struct
+{
+    volatile int lock;
+
+}Hash_table_spin_lock;
+
+typedef struct
+{
+	Count_Table** sub_h;
+    Hash_table_spin_lock* sub_h_lock;
+    int prefix_bits;
+    int suffix_bits;
+	int size;
+    uint64_t suffix_mode;
+} Total_Count_Table;
+
+/********************************for debug***************************************/
+inline void print_64bit(uint64_t x)
+{
+    int i;
+    for(i = 63; i >= 0; i--)
+	{
+		if(x & ((1ULL<<i)))
+			fprintf(stderr, "1");
+		else
+			fprintf(stderr, "0");
+	}
+
+    fprintf(stderr, "\n");
+}
+/********************************for debug***************************************/
+
+inline void get_sub_table(uint64_t* get_sub_ID, uint64_t* get_sub_key, Total_Count_Table* TCB, Hash_code* code, int k)
+{
+    uint64_t h_key, low_key;
+    ///k有可能是64，所以可能会有问题
+    ///low_key = code->x[0] | (code->x[1] << k);
+    low_key = code->x[0] | (code->x[1] << SAFE_SHIFT(k));
+    //k不可能为0, 所以这个右移不会有问题
+    h_key = code->x[1] >> (64 - k);
+
+    ///注意suffix_bits最大就是64
+    ///前一个右移不安全，因为TCB->suffix_bits有可能为64
+    ///后一个左移安全，因为TCB->suffix_bits不可能为0
+    //uint64_t sub_ID = (low_key >> TCB->suffix_bits) | (h_key << (64 - TCB->suffix_bits));
+    uint64_t sub_ID = (low_key >> SAFE_SHIFT(TCB->suffix_bits)) | (h_key << (64 - TCB->suffix_bits));
+    uint64_t sub_key = (low_key & TCB->suffix_mode);
+
+    *get_sub_ID = sub_ID;
+    *get_sub_key = sub_key;
+
+}
+
+
+inline void insert_Total_Count_Table(Total_Count_Table* TCB, Hash_code* code, int k)
+{
+    uint64_t sub_ID, sub_key;
+    get_sub_table(&sub_ID, &sub_key, TCB, code, k);
+
+    khint_t t;  ///这就是个迭代器
+    int absent;
+
+
+    while (__sync_lock_test_and_set(&TCB->sub_h_lock[sub_ID].lock, 1))
+    {
+        while (TCB->sub_h_lock[sub_ID].lock);
+    }
+
+    t = kh_put(COUNT64, TCB->sub_h[sub_ID], sub_key, &absent);
+    if (absent)
+    {
+        kh_value(TCB->sub_h[sub_ID], t) = 1;
+    }
+    else   ///哈希表中已有的元素
+    {
+        //kh_value(TCB->sub_h[sub_ID], t) = kh_value(TCB->sub_h[sub_ID], t) + 1;
+        kh_value(TCB->sub_h[sub_ID], t)++;
+    }
+
+    __sync_lock_release(&TCB->sub_h_lock[sub_ID].lock);
+}
+
+inline int get_Total_Count_Table(Total_Count_Table* TCB, Hash_code* code, int k)
+{
+
+    uint64_t sub_ID, sub_key;
+    get_sub_table(&sub_ID, &sub_key, TCB, code, k);
+
+    khint_t t;  ///这就是个迭代器
+    int absent;
+
+    ///查询哈希表，key为k
+    t = kh_get(COUNT64, TCB->sub_h[sub_ID], sub_key);
+
+    if (t != kh_end(TCB->sub_h[sub_ID]))
+    {
+        return kh_value(TCB->sub_h[sub_ID], t);
+    }
+    else
+    {
+        return -1;
+    }
+
+}
+
+/********************************for debug***************************************/
+inline int verify_Total_Count_Table(Total_Count_Table* TCB, Hash_code* code, int k)
+{
+    uint64_t sub_ID, sub_key;
+    get_sub_table(&sub_ID, &sub_key, TCB, code, k);
+
+    khint_t t;  ///这就是个迭代器
+    int absent;
+
+    ///查询哈希表，key为k
+    t = kh_get(COUNT64, TCB->sub_h[sub_ID], sub_key);
+
+    if (t != kh_end(TCB->sub_h[sub_ID]))
+    {
+        kh_value(TCB->sub_h[sub_ID], t)--;
+        if (kh_value(TCB->sub_h[sub_ID], t)<0)
+        {
+           return -1;
+        }
+        else
+        {
+            return 1;
+        }
+    }
+    else
+    {
+        return -1;
+    }
+}
+/********************************for debug***************************************/
+
+/********************************for debug***************************************/
+inline int Traverse_Total_Count_Table(Total_Count_Table* TCB)
+{
+    int i;
+    Count_Table* h;
+    khint_t k;
+
+    long long non_empty_k_mer = 0;
+
+    for (i = 0; i < TCB->size; i++)
+    {
+        h = TCB->sub_h[i];
+        for (k = kh_begin(h); k != kh_end(h); ++k)
+        {
+            if (kh_exist(h, k))            // test if a bucket contains data
+            {
+                non_empty_k_mer++;
+
+                if (kh_value(h, k)!= 0)
+                {
+                    fprintf(stderr, "ERROR when Traversing!\n");
+                }
+            }
+        }
+    }
+    
+    fprintf(stdout, "non_empty_k_mer: %lld\n", non_empty_k_mer);
+}
+/********************************for debug***************************************/
+
+
+void init_Total_Count_Table(int k, Total_Count_Table* TCB);
+void destory_Total_Count_Table(Total_Count_Table* TCB);
+
+void init_Count_Table(Count_Table** table);
+
+/********************************for debug***************************************/
+void test_COUNT64();
+/********************************for debug***************************************/
+
+#endif

Makefile

@@ -3,7 +3,7 @@ CC=g++
 CFLAGS = -w -c -msse4.2 -mpopcnt -fomit-frame-pointer -Winline -O3 -lz
 LDFLAGS = -lm -lz -lpthread -O3 -mpopcnt -msse4.2 -lz -w
 
-SOURCES = main.cpp CommandLines.cpp Process_Read.cpp Assembly.cpp
+SOURCES = main.cpp CommandLines.cpp Process_Read.cpp Assembly.cpp kmer.cpp Hash_Table.cpp
 OBJECTS = $(SOURCES:.c=.o)
 EXECUTABLE = ccs_assembly
 

Process_Read.cpp

@@ -65,13 +65,8 @@ int get_read(kseq_t *s)
 
 void init_R_buffer_block(R_buffer_block* curr_sub_block)
 {
-	int i;
-
-	for (i = 0; i < RDB.size; i++)
-	{
-		curr_sub_block->read = (kseq_t*)calloc(RDB.block_inner_size, sizeof(kseq_t));	
-		curr_sub_block->num = 0;
-	}	
+	curr_sub_block->read = (kseq_t*)calloc(RDB.block_inner_size, sizeof(kseq_t));	
+	curr_sub_block->num = 0;
 }
 
 
@@ -94,6 +89,27 @@ void init_R_buffer(int thread_num)
 }
 
 
+void destory_R_buffer_block(R_buffer_block* curr_sub_block)
+{
+
+	free(curr_sub_block->read);
+}
+
+
+void destory_R_buffer()
+{
+	int i = 0;
+
+	for (i = 0; i < RDB.size; i++)
+	{
+		destory_R_buffer_block(&RDB.sub_block[i]);
+	}
+
+	free(RDB.sub_block);
+	
+}
+
+
 inline void load_read_block(R_buffer_block* read_batch, int batch_read_size,
 	int* return_file_flag)
 {
@@ -215,6 +231,8 @@ void* input_reads_muti_threads(void*)
 	pthread_cond_signal(&i_readinputstallCond);  //important
 	pthread_mutex_unlock(&i_readinputMutex);
 
+	destory_R_buffer_block(&tmp_buf);
+
 }
 
 

Process_Read.h

@@ -8,7 +8,7 @@
 #include "kseq.h"
 
 #define READ_BLOCK_SIZE 64
-#define READ_BLOCK_NUM_PRE_THR 32
+#define READ_BLOCK_NUM_PRE_THR 100
 
 #define IS_FULL(buffer) ((buffer.num >= buffer.size)?1:0)
 #define IS_EMPTY(buffer) ((buffer.num == 0)?1:0)
@@ -46,5 +46,7 @@ void init_R_buffer_block(R_buffer_block* curr_sub_block);
 int get_reads_mul_thread(R_buffer_block* curr_sub_block);
 
 void Counting_block();
+void destory_R_buffer_block(R_buffer_block* curr_sub_block);
+void destory_R_buffer();
 
 #endif

khash.h

@@ -0,0 +1,627 @@
+/* The MIT License
+
+   Copyright (c) 2008, 2009, 2011 by Attractive Chaos <attractor@live.co.uk>
+
+   Permission is hereby granted, free of charge, to any person obtaining
+   a copy of this software and associated documentation files (the
+   "Software"), to deal in the Software without restriction, including
+   without limitation the rights to use, copy, modify, merge, publish,
+   distribute, sublicense, and/or sell copies of the Software, and to
+   permit persons to whom the Software is furnished to do so, subject to
+   the following conditions:
+
+   The above copyright notice and this permission notice shall be
+   included in all copies or substantial portions of the Software.
+
+   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+   EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+   MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+   NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
+   BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
+   ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
+   CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+   SOFTWARE.
+*/
+
+/*
+  An example:
+
+#include "khash.h"
+KHASH_MAP_INIT_INT(32, char)
+int main() {
+	int ret, is_missing;
+	khiter_t k;
+	khash_t(32) *h = kh_init(32);
+	k = kh_put(32, h, 5, &ret);
+	kh_value(h, k) = 10;
+	k = kh_get(32, h, 10);
+	is_missing = (k == kh_end(h));
+	k = kh_get(32, h, 5);
+	kh_del(32, h, k);
+	for (k = kh_begin(h); k != kh_end(h); ++k)
+		if (kh_exist(h, k)) kh_value(h, k) = 1;
+	kh_destroy(32, h);
+	return 0;
+}
+*/
+
+/*
+  2013-05-02 (0.2.8):
+
+	* Use quadratic probing. When the capacity is power of 2, stepping function
+	  i*(i+1)/2 guarantees to traverse each bucket. It is better than double
+	  hashing on cache performance and is more robust than linear probing.
+
+	  In theory, double hashing should be more robust than quadratic probing.
+	  However, my implementation is probably not for large hash tables, because
+	  the second hash function is closely tied to the first hash function,
+	  which reduce the effectiveness of double hashing.
+
+	Reference: http://research.cs.vt.edu/AVresearch/hashing/quadratic.php
+
+  2011-12-29 (0.2.7):
+
+    * Minor code clean up; no actual effect.
+
+  2011-09-16 (0.2.6):
+
+	* The capacity is a power of 2. This seems to dramatically improve the
+	  speed for simple keys. Thank Zilong Tan for the suggestion. Reference:
+
+	   - http://code.google.com/p/ulib/
+	   - http://nothings.org/computer/judy/
+
+	* Allow to optionally use linear probing which usually has better
+	  performance for random input. Double hashing is still the default as it
+	  is more robust to certain non-random input.
+
+	* Added Wang's integer hash function (not used by default). This hash
+	  function is more robust to certain non-random input.
+
+  2011-02-14 (0.2.5):
+
+    * Allow to declare global functions.
+
+  2009-09-26 (0.2.4):
+
+    * Improve portability
+
+  2008-09-19 (0.2.3):
+
+	* Corrected the example
+	* Improved interfaces
+
+  2008-09-11 (0.2.2):
+
+	* Improved speed a little in kh_put()
+
+  2008-09-10 (0.2.1):
+
+	* Added kh_clear()
+	* Fixed a compiling error
+
+  2008-09-02 (0.2.0):
+
+	* Changed to token concatenation which increases flexibility.
+
+  2008-08-31 (0.1.2):
+
+	* Fixed a bug in kh_get(), which has not been tested previously.
+
+  2008-08-31 (0.1.1):
+
+	* Added destructor
+*/
+
+
+#ifndef __AC_KHASH_H
+#define __AC_KHASH_H
+
+/*!
+  @header
+
+  Generic hash table library.
+ */
+
+#define AC_VERSION_KHASH_H "0.2.8"
+
+#include <stdlib.h>
+#include <string.h>
+#include <limits.h>
+
+/* compiler specific configuration */
+
+#if UINT_MAX == 0xffffffffu
+typedef unsigned int khint32_t;
+#elif ULONG_MAX == 0xffffffffu
+typedef unsigned long khint32_t;
+#endif
+
+#if ULONG_MAX == ULLONG_MAX
+typedef unsigned long khint64_t;
+#else
+typedef unsigned long long khint64_t;
+#endif
+
+#ifndef kh_inline
+#ifdef _MSC_VER
+#define kh_inline __inline
+#else
+#define kh_inline inline
+#endif
+#endif /* kh_inline */
+
+#ifndef klib_unused
+#if (defined __clang__ && __clang_major__ >= 3) || (defined __GNUC__ && __GNUC__ >= 3)
+#define klib_unused __attribute__ ((__unused__))
+#else
+#define klib_unused
+#endif
+#endif /* klib_unused */
+
+typedef khint32_t khint_t;
+typedef khint_t khiter_t;
+
+#define __ac_isempty(flag, i) ((flag[i>>4]>>((i&0xfU)<<1))&2)
+#define __ac_isdel(flag, i) ((flag[i>>4]>>((i&0xfU)<<1))&1)
+#define __ac_iseither(flag, i) ((flag[i>>4]>>((i&0xfU)<<1))&3)
+#define __ac_set_isdel_false(flag, i) (flag[i>>4]&=~(1ul<<((i&0xfU)<<1)))
+#define __ac_set_isempty_false(flag, i) (flag[i>>4]&=~(2ul<<((i&0xfU)<<1)))
+#define __ac_set_isboth_false(flag, i) (flag[i>>4]&=~(3ul<<((i&0xfU)<<1)))
+#define __ac_set_isdel_true(flag, i) (flag[i>>4]|=1ul<<((i&0xfU)<<1))
+
+#define __ac_fsize(m) ((m) < 16? 1 : (m)>>4)
+
+#ifndef kroundup32
+#define kroundup32(x) (--(x), (x)|=(x)>>1, (x)|=(x)>>2, (x)|=(x)>>4, (x)|=(x)>>8, (x)|=(x)>>16, ++(x))
+#endif
+
+#ifndef kcalloc
+#define kcalloc(N,Z) calloc(N,Z)
+#endif
+#ifndef kmalloc
+#define kmalloc(Z) malloc(Z)
+#endif
+#ifndef krealloc
+#define krealloc(P,Z) realloc(P,Z)
+#endif
+#ifndef kfree
+#define kfree(P) free(P)
+#endif
+
+static const double __ac_HASH_UPPER = 0.77;
+
+#define __KHASH_TYPE(name, khkey_t, khval_t) \
+	typedef struct kh_##name##_s { \
+		khint_t n_buckets, size, n_occupied, upper_bound; \
+		khint32_t *flags; \
+		khkey_t *keys; \
+		khval_t *vals; \
+	} kh_##name##_t;
+
+#define __KHASH_PROTOTYPES(name, khkey_t, khval_t)	 					\
+	extern kh_##name##_t *kh_init_##name(void);							\
+	extern void kh_destroy_##name(kh_##name##_t *h);					\
+	extern void kh_clear_##name(kh_##name##_t *h);						\
+	extern khint_t kh_get_##name(const kh_##name##_t *h, khkey_t key); 	\
+	extern int kh_resize_##name(kh_##name##_t *h, khint_t new_n_buckets); \
+	extern khint_t kh_put_##name(kh_##name##_t *h, khkey_t key, int *ret); \
+	extern void kh_del_##name(kh_##name##_t *h, khint_t x);
+
+#define __KHASH_IMPL(name, SCOPE, khkey_t, khval_t, kh_is_map, __hash_func, __hash_equal) \
+	SCOPE kh_##name##_t *kh_init_##name(void) {							\
+		return (kh_##name##_t*)kcalloc(1, sizeof(kh_##name##_t));		\
+	}																	\
+	SCOPE void kh_destroy_##name(kh_##name##_t *h)						\
+	{																	\
+		if (h) {														\
+			kfree((void *)h->keys); kfree(h->flags);					\
+			kfree((void *)h->vals);										\
+			kfree(h);													\
+		}																\
+	}																	\
+	SCOPE void kh_clear_##name(kh_##name##_t *h)						\
+	{																	\
+		if (h && h->flags) {											\
+			memset(h->flags, 0xaa, __ac_fsize(h->n_buckets) * sizeof(khint32_t)); \
+			h->size = h->n_occupied = 0;								\
+		}																\
+	}																	\
+	SCOPE khint_t kh_get_##name(const kh_##name##_t *h, khkey_t key) 	\
+	{																	\
+		if (h->n_buckets) {												\
+			khint_t k, i, last, mask, step = 0; \
+			mask = h->n_buckets - 1;									\
+			k = __hash_func(key); i = k & mask;							\
+			last = i; \
+			while (!__ac_isempty(h->flags, i) && (__ac_isdel(h->flags, i) || !__hash_equal(h->keys[i], key))) { \
+				i = (i + (++step)) & mask; \
+				if (i == last) return h->n_buckets;						\
+			}															\
+			return __ac_iseither(h->flags, i)? h->n_buckets : i;		\
+		} else return 0;												\
+	}																	\
+	SCOPE int kh_resize_##name(kh_##name##_t *h, khint_t new_n_buckets) \
+	{ /* This function uses 0.25*n_buckets bytes of working space instead of [sizeof(key_t+val_t)+.25]*n_buckets. */ \
+		khint32_t *new_flags = 0;										\
+		khint_t j = 1;													\
+		{																\
+			kroundup32(new_n_buckets); 									\
+			if (new_n_buckets < 4) new_n_buckets = 4;					\
+			if (h->size >= (khint_t)(new_n_buckets * __ac_HASH_UPPER + 0.5)) j = 0;	/* requested size is too small */ \
+			else { /* hash table size to be changed (shrink or expand); rehash */ \
+				new_flags = (khint32_t*)kmalloc(__ac_fsize(new_n_buckets) * sizeof(khint32_t));	\
+				if (!new_flags) return -1;								\
+				memset(new_flags, 0xaa, __ac_fsize(new_n_buckets) * sizeof(khint32_t)); \
+				if (h->n_buckets < new_n_buckets) {	/* expand */		\
+					khkey_t *new_keys = (khkey_t*)krealloc((void *)h->keys, new_n_buckets * sizeof(khkey_t)); \
+					if (!new_keys) { kfree(new_flags); return -1; }		\
+					h->keys = new_keys;									\
+					if (kh_is_map) {									\
+						khval_t *new_vals = (khval_t*)krealloc((void *)h->vals, new_n_buckets * sizeof(khval_t)); \
+						if (!new_vals) { kfree(new_flags); return -1; }	\
+						h->vals = new_vals;								\
+					}													\
+				} /* otherwise shrink */								\
+			}															\
+		}																\
+		if (j) { /* rehashing is needed */								\
+			for (j = 0; j != h->n_buckets; ++j) {						\
+				if (__ac_iseither(h->flags, j) == 0) {					\
+					khkey_t key = h->keys[j];							\
+					khval_t val;										\
+					khint_t new_mask;									\
+					new_mask = new_n_buckets - 1; 						\
+					if (kh_is_map) val = h->vals[j];					\
+					__ac_set_isdel_true(h->flags, j);					\
+					while (1) { /* kick-out process; sort of like in Cuckoo hashing */ \
+						khint_t k, i, step = 0; \
+						k = __hash_func(key);							\
+						i = k & new_mask;								\
+						while (!__ac_isempty(new_flags, i)) i = (i + (++step)) & new_mask; \
+						__ac_set_isempty_false(new_flags, i);			\
+						if (i < h->n_buckets && __ac_iseither(h->flags, i) == 0) { /* kick out the existing element */ \
+							{ khkey_t tmp = h->keys[i]; h->keys[i] = key; key = tmp; } \
+							if (kh_is_map) { khval_t tmp = h->vals[i]; h->vals[i] = val; val = tmp; } \
+							__ac_set_isdel_true(h->flags, i); /* mark it as deleted in the old hash table */ \
+						} else { /* write the element and jump out of the loop */ \
+							h->keys[i] = key;							\
+							if (kh_is_map) h->vals[i] = val;			\
+							break;										\
+						}												\
+					}													\
+				}														\
+			}															\
+			if (h->n_buckets > new_n_buckets) { /* shrink the hash table */ \
+				h->keys = (khkey_t*)krealloc((void *)h->keys, new_n_buckets * sizeof(khkey_t)); \
+				if (kh_is_map) h->vals = (khval_t*)krealloc((void *)h->vals, new_n_buckets * sizeof(khval_t)); \
+			}															\
+			kfree(h->flags); /* free the working space */				\
+			h->flags = new_flags;										\
+			h->n_buckets = new_n_buckets;								\
+			h->n_occupied = h->size;									\
+			h->upper_bound = (khint_t)(h->n_buckets * __ac_HASH_UPPER + 0.5); \
+		}																\
+		return 0;														\
+	}																	\
+	SCOPE khint_t kh_put_##name(kh_##name##_t *h, khkey_t key, int *ret) \
+	{																	\
+		khint_t x;														\
+		if (h->n_occupied >= h->upper_bound) { /* update the hash table */ \
+			if (h->n_buckets > (h->size<<1)) {							\
+				if (kh_resize_##name(h, h->n_buckets - 1) < 0) { /* clear "deleted" elements */ \
+					*ret = -1; return h->n_buckets;						\
+				}														\
+			} else if (kh_resize_##name(h, h->n_buckets + 1) < 0) { /* expand the hash table */ \
+				*ret = -1; return h->n_buckets;							\
+			}															\
+		} /* TODO: to implement automatically shrinking; resize() already support shrinking */ \
+		{																\
+			khint_t k, i, site, last, mask = h->n_buckets - 1, step = 0; \
+			x = site = h->n_buckets; k = __hash_func(key); i = k & mask; \
+			if (__ac_isempty(h->flags, i)) x = i; /* for speed up */	\
+			else {														\
+				last = i; \
+				while (!__ac_isempty(h->flags, i) && (__ac_isdel(h->flags, i) || !__hash_equal(h->keys[i], key))) { \
+					if (__ac_isdel(h->flags, i)) site = i;				\
+					i = (i + (++step)) & mask; \
+					if (i == last) { x = site; break; }					\
+				}														\
+				if (x == h->n_buckets) {								\
+					if (__ac_isempty(h->flags, i) && site != h->n_buckets) x = site; \
+					else x = i;											\
+				}														\
+			}															\
+		}																\
+		if (__ac_isempty(h->flags, x)) { /* not present at all */		\
+			h->keys[x] = key;											\
+			__ac_set_isboth_false(h->flags, x);							\
+			++h->size; ++h->n_occupied;									\
+			*ret = 1;													\
+		} else if (__ac_isdel(h->flags, x)) { /* deleted */				\
+			h->keys[x] = key;											\
+			__ac_set_isboth_false(h->flags, x);							\
+			++h->size;													\
+			*ret = 2;													\
+		} else *ret = 0; /* Don't touch h->keys[x] if present and not deleted */ \
+		return x;														\
+	}																	\
+	SCOPE void kh_del_##name(kh_##name##_t *h, khint_t x)				\
+	{																	\
+		if (x != h->n_buckets && !__ac_iseither(h->flags, x)) {			\
+			__ac_set_isdel_true(h->flags, x);							\
+			--h->size;													\
+		}																\
+	}
+
+#define KHASH_DECLARE(name, khkey_t, khval_t)		 					\
+	__KHASH_TYPE(name, khkey_t, khval_t) 								\
+	__KHASH_PROTOTYPES(name, khkey_t, khval_t)
+
+#define KHASH_INIT2(name, SCOPE, khkey_t, khval_t, kh_is_map, __hash_func, __hash_equal) \
+	__KHASH_TYPE(name, khkey_t, khval_t) 								\
+	__KHASH_IMPL(name, SCOPE, khkey_t, khval_t, kh_is_map, __hash_func, __hash_equal)
+
+#define KHASH_INIT(name, khkey_t, khval_t, kh_is_map, __hash_func, __hash_equal) \
+	KHASH_INIT2(name, static kh_inline klib_unused, khkey_t, khval_t, kh_is_map, __hash_func, __hash_equal)
+
+/* --- BEGIN OF HASH FUNCTIONS --- */
+
+/*! @function
+  @abstract     Integer hash function
+  @param  key   The integer [khint32_t]
+  @return       The hash value [khint_t]
+ */
+#define kh_int_hash_func(key) (khint32_t)(key)
+/*! @function
+  @abstract     Integer comparison function
+ */
+#define kh_int_hash_equal(a, b) ((a) == (b))
+/*! @function
+  @abstract     64-bit integer hash function
+  @param  key   The integer [khint64_t]
+  @return       The hash value [khint_t]
+ */
+#define kh_int64_hash_func(key) (khint32_t)((key)>>33^(key)^(key)<<11)
+/*! @function
+  @abstract     64-bit integer comparison function
+ */
+#define kh_int64_hash_equal(a, b) ((a) == (b))
+/*! @function
+  @abstract     const char* hash function
+  @param  s     Pointer to a null terminated string
+  @return       The hash value
+ */
+static kh_inline khint_t __ac_X31_hash_string(const char *s)
+{
+	khint_t h = (khint_t)*s;
+	if (h) for (++s ; *s; ++s) h = (h << 5) - h + (khint_t)*s;
+	return h;
+}
+/*! @function
+  @abstract     Another interface to const char* hash function
+  @param  key   Pointer to a null terminated string [const char*]
+  @return       The hash value [khint_t]
+ */
+#define kh_str_hash_func(key) __ac_X31_hash_string(key)
+/*! @function
+  @abstract     Const char* comparison function
+ */
+#define kh_str_hash_equal(a, b) (strcmp(a, b) == 0)
+
+static kh_inline khint_t __ac_Wang_hash(khint_t key)
+{
+    key += ~(key << 15);
+    key ^=  (key >> 10);
+    key +=  (key << 3);
+    key ^=  (key >> 6);
+    key += ~(key << 11);
+    key ^=  (key >> 16);
+    return key;
+}
+#define kh_int_hash_func2(key) __ac_Wang_hash((khint_t)key)
+
+/* --- END OF HASH FUNCTIONS --- */
+
+/* Other convenient macros... */
+
+/*!
+  @abstract Type of the hash table.
+  @param  name  Name of the hash table [symbol]
+ */
+#define khash_t(name) kh_##name##_t
+
+/*! @function
+  @abstract     Initiate a hash table.
+  @param  name  Name of the hash table [symbol]
+  @return       Pointer to the hash table [khash_t(name)*]
+ */
+#define kh_init(name) kh_init_##name()
+
+/*! @function
+  @abstract     Destroy a hash table.
+  @param  name  Name of the hash table [symbol]
+  @param  h     Pointer to the hash table [khash_t(name)*]
+ */
+#define kh_destroy(name, h) kh_destroy_##name(h)
+
+/*! @function
+  @abstract     Reset a hash table without deallocating memory.
+  @param  name  Name of the hash table [symbol]
+  @param  h     Pointer to the hash table [khash_t(name)*]
+ */
+#define kh_clear(name, h) kh_clear_##name(h)
+
+/*! @function
+  @abstract     Resize a hash table.
+  @param  name  Name of the hash table [symbol]
+  @param  h     Pointer to the hash table [khash_t(name)*]
+  @param  s     New size [khint_t]
+ */
+#define kh_resize(name, h, s) kh_resize_##name(h, s)
+
+/*! @function
+  @abstract     Insert a key to the hash table.
+  @param  name  Name of the hash table [symbol]
+  @param  h     Pointer to the hash table [khash_t(name)*]
+  @param  k     Key [type of keys]
+  @param  r     Extra return code: -1 if the operation failed;
+                0 if the key is present in the hash table;
+                1 if the bucket is empty (never used); 2 if the element in
+				the bucket has been deleted [int*]
+  @return       Iterator to the inserted element [khint_t]
+ */
+#define kh_put(name, h, k, r) kh_put_##name(h, k, r)
+
+/*! @function
+  @abstract     Retrieve a key from the hash table.
+  @param  name  Name of the hash table [symbol]
+  @param  h     Pointer to the hash table [khash_t(name)*]
+  @param  k     Key [type of keys]
+  @return       Iterator to the found element, or kh_end(h) if the element is absent [khint_t]
+ */
+#define kh_get(name, h, k) kh_get_##name(h, k)
+
+/*! @function
+  @abstract     Remove a key from the hash table.
+  @param  name  Name of the hash table [symbol]
+  @param  h     Pointer to the hash table [khash_t(name)*]
+  @param  k     Iterator to the element to be deleted [khint_t]
+ */
+#define kh_del(name, h, k) kh_del_##name(h, k)
+
+/*! @function
+  @abstract     Test whether a bucket contains data.
+  @param  h     Pointer to the hash table [khash_t(name)*]
+  @param  x     Iterator to the bucket [khint_t]
+  @return       1 if containing data; 0 otherwise [int]
+ */
+#define kh_exist(h, x) (!__ac_iseither((h)->flags, (x)))
+
+/*! @function
+  @abstract     Get key given an iterator
+  @param  h     Pointer to the hash table [khash_t(name)*]
+  @param  x     Iterator to the bucket [khint_t]
+  @return       Key [type of keys]
+ */
+#define kh_key(h, x) ((h)->keys[x])
+
+/*! @function
+  @abstract     Get value given an iterator
+  @param  h     Pointer to the hash table [khash_t(name)*]
+  @param  x     Iterator to the bucket [khint_t]
+  @return       Value [type of values]
+  @discussion   For hash sets, calling this results in segfault.
+ */
+#define kh_val(h, x) ((h)->vals[x])
+
+/*! @function
+  @abstract     Alias of kh_val()
+ */
+#define kh_value(h, x) ((h)->vals[x])
+
+/*! @function
+  @abstract     Get the start iterator
+  @param  h     Pointer to the hash table [khash_t(name)*]
+  @return       The start iterator [khint_t]
+ */
+#define kh_begin(h) (khint_t)(0)
+
+/*! @function
+  @abstract     Get the end iterator
+  @param  h     Pointer to the hash table [khash_t(name)*]
+  @return       The end iterator [khint_t]
+ */
+#define kh_end(h) ((h)->n_buckets)
+
+/*! @function
+  @abstract     Get the number of elements in the hash table
+  @param  h     Pointer to the hash table [khash_t(name)*]
+  @return       Number of elements in the hash table [khint_t]
+ */
+#define kh_size(h) ((h)->size)
+
+/*! @function
+  @abstract     Get the number of buckets in the hash table
+  @param  h     Pointer to the hash table [khash_t(name)*]
+  @return       Number of buckets in the hash table [khint_t]
+ */
+#define kh_n_buckets(h) ((h)->n_buckets)
+
+/*! @function
+  @abstract     Iterate over the entries in the hash table
+  @param  h     Pointer to the hash table [khash_t(name)*]
+  @param  kvar  Variable to which key will be assigned
+  @param  vvar  Variable to which value will be assigned
+  @param  code  Block of code to execute
+ */
+#define kh_foreach(h, kvar, vvar, code) { khint_t __i;		\
+	for (__i = kh_begin(h); __i != kh_end(h); ++__i) {		\
+		if (!kh_exist(h,__i)) continue;						\
+		(kvar) = kh_key(h,__i);								\
+		(vvar) = kh_val(h,__i);								\
+		code;												\
+	} }
+
+/*! @function
+  @abstract     Iterate over the values in the hash table
+  @param  h     Pointer to the hash table [khash_t(name)*]
+  @param  vvar  Variable to which value will be assigned
+  @param  code  Block of code to execute
+ */
+#define kh_foreach_value(h, vvar, code) { khint_t __i;		\
+	for (__i = kh_begin(h); __i != kh_end(h); ++__i) {		\
+		if (!kh_exist(h,__i)) continue;						\
+		(vvar) = kh_val(h,__i);								\
+		code;												\
+	} }
+
+/* More convenient interfaces */
+
+/*! @function
+  @abstract     Instantiate a hash set containing integer keys
+  @param  name  Name of the hash table [symbol]
+ */
+#define KHASH_SET_INIT_INT(name)										\
+	KHASH_INIT(name, khint32_t, char, 0, kh_int_hash_func, kh_int_hash_equal)
+
+/*! @function
+  @abstract     Instantiate a hash map containing integer keys
+  @param  name  Name of the hash table [symbol]
+  @param  khval_t  Type of values [type]
+ */
+#define KHASH_MAP_INIT_INT(name, khval_t)								\
+	KHASH_INIT(name, khint32_t, khval_t, 1, kh_int_hash_func, kh_int_hash_equal)
+
+/*! @function
+  @abstract     Instantiate a hash set containing 64-bit integer keys
+  @param  name  Name of the hash table [symbol]
+ */
+#define KHASH_SET_INIT_INT64(name)										\
+	KHASH_INIT(name, khint64_t, char, 0, kh_int64_hash_func, kh_int64_hash_equal)
+
+/*! @function
+  @abstract     Instantiate a hash map containing 64-bit integer keys
+  @param  name  Name of the hash table [symbol]
+  @param  khval_t  Type of values [type]
+ */
+#define KHASH_MAP_INIT_INT64(name, khval_t)								\
+	KHASH_INIT(name, khint64_t, khval_t, 1, kh_int64_hash_func, kh_int64_hash_equal)
+
+typedef const char *kh_cstr_t;
+/*! @function
+  @abstract     Instantiate a hash map containing const char* keys
+  @param  name  Name of the hash table [symbol]
+ */
+#define KHASH_SET_INIT_STR(name)										\
+	KHASH_INIT(name, kh_cstr_t, char, 0, kh_str_hash_func, kh_str_hash_equal)
+
+/*! @function
+  @abstract     Instantiate a hash map containing const char* keys
+  @param  name  Name of the hash table [symbol]
+  @param  khval_t  Type of values [type]
+ */
+#define KHASH_MAP_INIT_STR(name, khval_t)								\
+	KHASH_INIT(name, kh_cstr_t, khval_t, 1, kh_str_hash_func, kh_str_hash_equal)
+
+#endif /* __AC_KHASH_H */

kmer.cpp

@@ -0,0 +1,16 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include "kmer.h"
+
+void init_HPC_seq(HPC_seq* seq, char* str, long long l)
+{
+    seq->i = 0;
+    seq->l = l;
+    seq->str = str;
+}
+
+void init_Hash_code(Hash_code* code)
+{
+    code->x[0] = 0;
+    code->x[1] = 0;
+}

kmer.h

@@ -0,0 +1,89 @@
+#ifndef __KMER__
+#define __KMER__
+#include "Process_Read.h"
+
+#define ALL (0xffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff)
+
+#define SAFE_SHIFT(k) k & ((k < 64)?ALL:0)
+
+
+static unsigned char seq_nt6_table[256] = {
+    5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,
+    5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,
+    5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,
+    5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,
+    5, 0, 5, 1,  5, 5, 5, 2,  5, 5, 5, 5,  5, 5, 5, 5,
+    5, 5, 5, 5,  3, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,
+    5, 0, 5, 1,  5, 5, 5, 2,  5, 5, 5, 5,  5, 5, 5, 5,
+    5, 5, 5, 5,  3, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,
+    5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,
+    5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,
+    5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,
+    5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,
+    5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,
+    5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,
+    5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,
+    5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5
+};
+
+
+
+typedef struct
+{
+    ///最大64-mer
+    ///x[0]低位
+    ///x[1]高位
+    uint64_t x[2];
+
+} Hash_code;
+
+
+
+typedef struct
+{
+	char* str;
+	long long l;
+    long long i;
+
+} HPC_seq;
+
+
+inline uint64_t get_HPC_code(HPC_seq* seq)
+{
+
+    if(seq->i < seq ->l)
+    {
+        char code = seq->str[seq->i];
+
+        for (; seq->i < seq->l; seq->i++)
+        {
+            if (seq->str[seq->i] != code)
+            {
+                break;
+            }
+        }
+
+        return (uint64_t)seq_nt6_table[(uint8_t)code];
+    }
+    else
+    {
+        ///end
+        return 6;
+    }
+    
+}
+
+inline void k_mer_append(Hash_code* code, uint64_t c, int k)
+{
+
+    uint64_t mask = ALL >> (64 -k);
+
+    code->x[0] = ((code->x[0]<<1) | (c&1))  & mask;
+    code->x[1] = ((code->x[1]<<1) | (c>>1)) & mask;
+}
+
+void init_HPC_seq(HPC_seq* seq, char* str, long long l);
+void init_Hash_code(Hash_code* code);
+
+
+#endif

main.cpp

@@ -4,6 +4,17 @@
 #include "Process_Read.h"
 #include "Assembly.h"
 
+/********************************for debug***************************************/
+///使用这个函数的时候，必须把Counting_multiple_thr()里的destory_Total_Count_Table(&TCB)注释掉
+void debug_Counting()
+{
+    init_kseq(read_file_name);
+    Verify_Counting();
+    fprintf(stderr, "debug over!\n");
+    destory_kseq();
+}
+/********************************for debug***************************************/
+
 
 int main(int argc, char *argv[])
 {
@@ -13,15 +24,15 @@ int main(int argc, char *argv[])
 
     init_kseq(read_file_name);
 
-
+    fprintf(stdout, "k-mer length: %d\n",k_mer_length);
 
     ///Counting();
     Counting_multiple_thr();
 
-
-
     destory_kseq();
 
+    ///debug_Counting();
+
 
     return 1;
 }