renamed files for code cleanup

Assembly.cpp

@@ -9,7 +9,7 @@
 #include "POA.h"
 #include "Correct.h"
 #include "Output.h"
-#include "yak.h"
+#include "htab.h"
 
 Total_Count_Table TCB;
 Total_Pos_Table PCB;

Hash_Table.h

@@ -2,7 +2,7 @@
 #define __HASHTABLE__
 #include "khashl.h"
 #include "kmer.h"
-#include "yak.h"
+#include "htab.h"
 
 KHASHL_MAP_INIT(static inline, Count_Table, ha_ct, uint64_t, int, kh_hash_dummy, kh_eq_generic)
 KHASHL_MAP_INIT(static inline, Pos_Table, ha_pt, uint64_t, uint64_t, kh_hash_dummy, kh_eq_generic)

Makefile

@@ -4,7 +4,7 @@ CPPFLAGS=
 INCLUDES=
 OBJS=		Output.o CommandLines.o Process_Read.o Assembly.o Hash_Table.o \
 			POA.o Correct.o Levenshtein_distance.o Overlaps.o Trio.o kthread.o \
-			yak-bbf.o yak-count.o hist.o sketch.o yak-sys.o
+			htab.o hist.o sketch.o sys.o
 EXE=		hifiasm
 LIBS=		-lz -lpthread -lm
 
@@ -33,30 +33,29 @@ depend:
 # DO NOT DELETE
 
 Assembly.o: Assembly.h CommandLines.h Process_Read.h kseq.h Overlaps.h kvec.h
-Assembly.o: kdq.h kmer.h Hash_Table.h khashl.h yak.h POA.h Correct.h
+Assembly.o: kdq.h kmer.h Hash_Table.h khashl.h htab.h POA.h Correct.h
 Assembly.o: Levenshtein_distance.h Output.h
 CommandLines.o: CommandLines.h ketopt.h
 Correct.o: Correct.h Hash_Table.h khashl.h kmer.h Process_Read.h kseq.h
-Correct.o: Overlaps.h kvec.h kdq.h CommandLines.h yak.h
+Correct.o: Overlaps.h kvec.h kdq.h CommandLines.h htab.h
 Correct.o: Levenshtein_distance.h POA.h Assembly.h
 Hash_Table.o: Hash_Table.h khashl.h kmer.h Process_Read.h kseq.h Overlaps.h
-Hash_Table.o: kvec.h kdq.h CommandLines.h yak.h Correct.h
+Hash_Table.o: kvec.h kdq.h CommandLines.h htab.h Correct.h
 Hash_Table.o: Levenshtein_distance.h POA.h ksort.h
 Levenshtein_distance.o: Levenshtein_distance.h
 Output.o: Output.h CommandLines.h
 Overlaps.o: Overlaps.h kvec.h kdq.h ksort.h Process_Read.h kseq.h
-Overlaps.o: CommandLines.h Hash_Table.h khashl.h kmer.h yak.h Correct.h
+Overlaps.o: CommandLines.h Hash_Table.h khashl.h kmer.h htab.h Correct.h
 Overlaps.o: Levenshtein_distance.h POA.h
 POA.o: POA.h Hash_Table.h khashl.h kmer.h Process_Read.h kseq.h Overlaps.h
-POA.o: kvec.h kdq.h CommandLines.h yak.h Correct.h Levenshtein_distance.h
+POA.o: kvec.h kdq.h CommandLines.h htab.h Correct.h Levenshtein_distance.h
 Process_Read.o: Process_Read.h kseq.h Overlaps.h kvec.h kdq.h CommandLines.h
 Trio.o: khashl.h kthread.h Process_Read.h kseq.h Overlaps.h kvec.h kdq.h
-Trio.o: CommandLines.h yak.h
+Trio.o: CommandLines.h htab.h
-hist.o: yak.h
+hist.o: htab.h
+htab.o: kthread.h khashl.h kseq.h htab.h CommandLines.h
 kthread.o: kthread.h
 main.o: CommandLines.h Process_Read.h kseq.h Overlaps.h kvec.h kdq.h
-main.o: Assembly.h Levenshtein_distance.h yak.h
+main.o: Assembly.h Levenshtein_distance.h htab.h
-sketch.o: kvec.h yak.h
+sketch.o: kvec.h htab.h
-yak-bbf.o: yak.h
+sys.o: htab.h
-yak-count.o: CommandLines.h yak.h khashl.h kthread.h kseq.h
-yak-sys.o: yak.h

Trio.cpp

@@ -6,7 +6,7 @@
 #include "khashl.h" // hash table
 #include "kthread.h"
 #include "Process_Read.h"
-#include "yak.h"
+#include "htab.h"
 #include "CommandLines.h"
 
 #define YAK_MAX_KMER     31

hist.cpp

@@ -1,5 +1,5 @@
 #include <stdio.h>
-#include "yak.h"
+#include "htab.h"
 
 static void yak_hist_line(int c, int x, int exceed, int64_t cnt)
 {

htab.cpp

@@ -7,7 +7,7 @@
 #include "kthread.h"
 #include "khashl.h"
 #include "kseq.h"
-#include "yak.h"
+#include "htab.h"
 #include "CommandLines.h"
 
 #define YAK_COUNTER_BITS 12
@@ -57,6 +57,54 @@ void yak_copt_init(yak_copt_t *o)
 	o->chunk_size = 10000000;
 }
 
+/************************
+ * Blocked bloom filter *
+ ************************/
+
+typedef struct {
+	int n_shift, n_hashes;
+	uint8_t *b;
+} yak_bf_t;
+
+yak_bf_t *yak_bf_init(int n_shift, int n_hashes)
+{
+	yak_bf_t *b;
+	void *ptr = 0;
+	if (n_shift + YAK_BLK_SHIFT > 64 || n_shift < YAK_BLK_SHIFT) return 0;
+	CALLOC(b, 1);
+	b->n_shift = n_shift;
+	b->n_hashes = n_hashes;
+	posix_memalign(&ptr, 1<<(YAK_BLK_SHIFT-3), 1ULL<<(n_shift-3));
+	b->b = (uint8_t*)ptr;
+	bzero(b->b, 1ULL<<(n_shift-3));
+	return b;
+}
+
+void yak_bf_destroy(yak_bf_t *b)
+{
+	if (b == 0) return;
+	free(b->b); free(b);
+}
+
+int yak_bf_insert(yak_bf_t *b, uint64_t hash)
+{
+	int x = b->n_shift - YAK_BLK_SHIFT;
+	uint64_t y = hash & ((1ULL<<x) - 1);
+	int h1 = hash >> x & YAK_BLK_MASK;
+	int h2 = hash >> b->n_shift & YAK_BLK_MASK;
+	uint8_t *p = &b->b[y<<(YAK_BLK_SHIFT-3)];
+	int i, z = h1, cnt = 0;
+	if ((h2&31) == 0) h2 = (h2 + 1) & YAK_BLK_MASK; // otherwise we may repeatedly use a few bits
+	for (i = 0; i < b->n_hashes; z = (z + h2) & YAK_BLK_MASK) {
+		uint8_t *q = &p[z>>3], u;
+		u = 1<<(z&7);
+		cnt += !!(*q & u);
+		*q |= u;
+		++i;
+	}
+	return cnt;
+}
+
 /********************
  * Count hash table *
  ********************/
@@ -306,7 +354,7 @@ typedef struct {
 	ha_mz1_t *b;
 } ch_buf_t;
 
-static inline void ch_insert_buf(ch_buf_t *buf, int p, uint64_t y) // insert a k-mer $y to a linear buffer
+static inline void ct_insert_buf(ch_buf_t *buf, int p, uint64_t y) // insert a k-mer $y to a linear buffer
 {
 	int pre = y & ((1<<p) - 1);
 	ch_buf_t *b = &buf[pre];
@@ -340,7 +388,7 @@ static void count_seq_buf(ch_buf_t *buf, int k, int p, int len, const char *seq)
 			x[2] = x[2] >> 1 | (uint64_t)(1 - (c&1))  << shift;
 			x[3] = x[3] >> 1 | (uint64_t)(1 - (c>>1)) << shift;
 			if (++l >= k)
-				ch_insert_buf(buf, p, yak_hash_long(x));
+				ct_insert_buf(buf, p, yak_hash_long(x));
 		} else l = 0, x[0] = x[1] = x[2] = x[3] = 0; // if there is an "N", restart
 	}
 }
@@ -358,7 +406,7 @@ static void count_seq_buf_HPC(ch_buf_t *buf, int k, int p, int len, const char *
 				x[2] = x[2] >> 1 | (uint64_t)(1 - (c&1))  << shift;
 				x[3] = x[3] >> 1 | (uint64_t)(1 - (c>>1)) << shift;
 				if (++l >= k)
-					ch_insert_buf(buf, p, yak_hash_long(x));
+					ct_insert_buf(buf, p, yak_hash_long(x));
 				last = c;
 			}
 		} else l = 0, last = -1, x[0] = x[1] = x[2] = x[3] = 0; // if there is an "N", restart
@@ -375,6 +423,8 @@ typedef struct { // global data structure for kt_pipeline()
 	const yak_copt_t *opt;
 	const void *flt_tab;
 	int create_new, is_store;
+	uint64_t batch_offset;
+	uint64_t n_base;
 	kseq_t *ks;
 	ha_ct_t *ct;
 	ha_pt_t *pt;
@@ -405,7 +455,7 @@ static void worker_for_mz(void *data, long i, int tid)
 	st_data_t *s = (st_data_t*)data;
 	ha_mz1_v *b = &s->mz_buf[tid];
 	s->mz_buf[tid].n = 0;
-	ha_sketch(s->seq[i], s->len[i], s->p->opt->w, s->p->opt->k, 0, s->p->opt->is_HPC, b, s->p->flt_tab);
+	ha_sketch(s->seq[i], s->len[i], s->p->opt->w, s->p->opt->k, s->p->batch_offset + i, s->p->opt->is_HPC, b, s->p->flt_tab);
 	s->mz[i].n = s->mz[i].m = b->n;
 	MALLOC(s->mz[i].a, b->n);
 	memcpy(s->mz[i].a, b->a, b->n * sizeof(ha_mz1_t));
@@ -421,7 +471,11 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip
 		s->p = p;
 		while ((ret = kseq_read(p->ks)) >= 0) {
 			int l = p->ks->seq.l;
-			if (l < p->opt->k) continue;
+			if (p->batch_offset + s->n_seq >= (1<<28) - 1) {
+				fprintf(stderr, "ERROR: this implementation supports no more than %d reads\n", (1<<28) - 1);
+				exit(1);
+			}
+			p->n_base += l;
 			if (s->n_seq == s->m_seq) {
 				s->m_seq = s->m_seq < 16? 16 : s->m_seq + (s->m_seq>>1);
 				REALLOC(s->len, s->m_seq);
@@ -473,7 +527,7 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip
 			} else {
 				for (i = 0; i < s->n_seq; ++i)
 					for (j = 0; j < s->mz[i].n; ++j)
-						ch_insert_buf(s->buf, p->opt->pre, s->mz[i].a[j].x);
+						ct_insert_buf(s->buf, p->opt->pre, s->mz[i].a[j].x);
 			}
 			for (i = 0; i < s->n_seq; ++i) {
 				free(s->mz[i].a);
@@ -494,6 +548,7 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip
 			free(s->buf[i].a);
 		}
 		p->ct->tot += n_ins;
+		p->batch_offset += s->n_seq;
 		free(s->buf);
 		fprintf(stderr, "[M::%s::%.3f*%.2f] processed %d sequences; %ld distinct k-mers in the hash table\n", __func__,
 				yak_realtime(), yak_cputime() / yak_realtime(), s->n_seq, (long)p->ct->tot);

htab.h

@@ -1,5 +1,5 @@
-#ifndef __YAK_H__
+#ifndef __HA_HTAB_H__
-#define __YAK_H__
+#define __HA_HTAB_H__
 #define __STDC_LIMIT_MACROS
 #include <stdint.h>
 
@@ -17,11 +17,6 @@ typedef struct {
 
 typedef struct { uint32_t n, m; ha_mz1_t *a; } ha_mz1_v;
 
-typedef struct {
-	int n_shift, n_hashes;
-	uint8_t *b;
-} yak_bf_t;
-
 extern const unsigned char seq_nt4_table[256];
 
 int ha_ft_isflt(const void *hh, uint64_t y);
@@ -37,10 +32,6 @@ long yak_peakrss(void);
 void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf);
 int yak_analyze_count(int n_cnt, const int64_t *cnt, int *peak_het);
 
-yak_bf_t *yak_bf_init(int n_shift, int n_hashes);
-void yak_bf_destroy(yak_bf_t *b);
-int yak_bf_insert(yak_bf_t *b, uint64_t hash);
-
 static inline uint64_t yak_hash64(uint64_t key, uint64_t mask) // invertible integer hash function
 {
 	key = (~key + (key << 21)) & mask; // key = (key << 21) - key - 1;

main.cpp

@@ -4,7 +4,7 @@
 #include "Process_Read.h"
 #include "Assembly.h"
 #include "Levenshtein_distance.h"
-#include "yak.h"
+#include "htab.h"
 
 int main(int argc, char *argv[])
 {

sketch.cpp

@@ -3,7 +3,7 @@
 #include <assert.h>
 #include <string.h>
 #include "kvec.h"
-#include "yak.h"
+#include "htab.h"
 
 typedef struct { // a simplified version of kdq
 	int front, count;

sys.cpp

@@ -1,6 +1,6 @@
 #include <sys/resource.h>
 #include <sys/time.h>
-#include "yak.h"
+#include "htab.h"
 
 int yak_verbose = 3;
 

yak-bbf.cpp

@@ -1,42 +0,0 @@
-#include <stdlib.h>
-#include <string.h>
-#include "yak.h"
-
-yak_bf_t *yak_bf_init(int n_shift, int n_hashes)
-{
-	yak_bf_t *b;
-	void *ptr = 0;
-	if (n_shift + YAK_BLK_SHIFT > 64 || n_shift < YAK_BLK_SHIFT) return 0;
-	CALLOC(b, 1);
-	b->n_shift = n_shift;
-	b->n_hashes = n_hashes;
-	posix_memalign(&ptr, 1<<(YAK_BLK_SHIFT-3), 1ULL<<(n_shift-3));
-	b->b = (uint8_t*)ptr;
-	bzero(b->b, 1ULL<<(n_shift-3));
-	return b;
-}
-
-void yak_bf_destroy(yak_bf_t *b)
-{
-	if (b == 0) return;
-	free(b->b); free(b);
-}
-
-int yak_bf_insert(yak_bf_t *b, uint64_t hash)
-{
-	int x = b->n_shift - YAK_BLK_SHIFT;
-	uint64_t y = hash & ((1ULL<<x) - 1);
-	int h1 = hash >> x & YAK_BLK_MASK;
-	int h2 = hash >> b->n_shift & YAK_BLK_MASK;
-	uint8_t *p = &b->b[y<<(YAK_BLK_SHIFT-3)];
-	int i, z = h1, cnt = 0;
-	if ((h2&31) == 0) h2 = (h2 + 1) & YAK_BLK_MASK; // otherwise we may repeatedly use a few bits
-	for (i = 0; i < b->n_hashes; z = (z + h2) & YAK_BLK_MASK) {
-		uint8_t *q = &p[z>>3], u;
-		u = 1<<(z&7);
-		cnt += !!(*q & u);
-		*q |= u;
-		++i;
-	}
-	return cnt;
-}