r192: estimate buffer memory

Assembly.cpp

@@ -429,11 +429,49 @@ void ha_ovec_destroy(ha_ovec_buf_t *b)
 	free(b);
 }
 
+static int64_t ha_Graph_mem(const Graph *g)
+{
+	int64_t i, mem = 0;
+	mem = sizeof(Graph) + g->node_q.size * 8 + g->g_nodes.size * sizeof(Node);
+	for (i = 0; i < (int64_t)g->g_nodes.size; ++i) {
+		Node *n = &g->g_nodes.list[i];
+		mem += n->mismatch_edges.size * sizeof(Edge);
+		mem += n->deletion_edges.size * sizeof(Edge);
+		mem += n->insertion_edges.size * sizeof(Edge);
+	}
+	mem += g->g_nodes.sort.size * 9;
+	return mem;
+}
+
+int64_t ha_ovec_mem(const ha_ovec_buf_t *b)
+{
+	int64_t i, mem = 0, mem_clist, mem_olist;
+	mem_clist = b->clist.size * sizeof(k_mer_hit) + b->clist.chainDP.size * 7 * 4;
+	mem_olist = b->olist.size * sizeof(overlap_region);
+	for (i = 0; i < (int64_t)b->olist.size; ++i) {
+		const overlap_region *r = &b->olist.list[i];
+		mem_olist += r->w_list_size * sizeof(window_list);
+		mem_olist += r->f_cigar.size * 8;
+		mem_olist += r->boundary_cigars.size * sizeof(window_list);
+	}
+	mem = ha_abuf_mem(b->ab) + mem_clist + mem_olist;
+	if (!b->is_final) {
+		mem += sizeof(Cigar_record) + b->cigar1.lost_base_size + b->cigar1.size * 4;
+		mem += sizeof(Correct_dumy) + b->correct.size * 8;
+		mem += sizeof(Round2_alignment) + b->round2.cigar.size * 4 + b->round2.tmp_cigar.size * 4;
+		mem += sizeof(haplotype_evdience_alloc) + b->hap.size * sizeof(haplotype_evdience) + b->hap.snp_matrix_size + b->hap.snp_stat_size * sizeof(SnpStats);
+		mem += ha_Graph_mem(&b->POA_Graph);
+		mem += ha_Graph_mem(&b->DAGCon);
+	}
+	return mem;
+}
+
 void* Overlap_calculate_heap_merge(void* arg)
 {
 	long long num_read_base = 0;
 	long long num_correct_base = 0;
 	long long num_recorrect_base = 0;
+	long long mem_buf;
 	int fully_cov, abnormal;
 
 	int thr_ID = *((int*)arg);
@@ -469,12 +507,14 @@ void* Overlap_calculate_heap_merge(void* arg)
 		push_overlaps(&(R_INF.reverse_paf[i]), &b->olist, 2, &R_INF, asm_opt.roundID%2);
 	}
 	finish_output_buffer();
+	mem_buf = ha_ovec_mem(b);
 	ha_ovec_destroy(b);
 
 	pthread_mutex_lock(&statistics);
 	asm_opt.num_bases += num_read_base;
 	asm_opt.num_corrected_bases += num_correct_base;
 	asm_opt.num_recorrected_bases += num_recorrect_base;
+	asm_opt.mem_buf += mem_buf;
 	pthread_mutex_unlock(&statistics);
 	return NULL;
 }
@@ -1179,7 +1219,7 @@ int ha_assemble(void)
 	int r, ovlp_loaded = 0;
 	if (asm_opt.load_index_from_disk && load_all_data_from_disk(&R_INF.paf, &R_INF.reverse_paf, asm_opt.output_file_name)) {
 		ovlp_loaded = 1;
-		fprintf(stderr, "[M::%s::%.3f*%.2f] ==> loaded overlaps from disk\n", __func__, yak_realtime(), yak_cputime() / yak_realtime());
+		fprintf(stderr, "[M::%s::%.3f*%.2f] ==> loaded overlaps from disk\n", __func__, yak_realtime(), yak_cpu_usage());
 	}
 	if (!ovlp_loaded) {
 		// construct hash table for high occurrence k-mers
@@ -1190,16 +1230,19 @@ int ha_assemble(void)
 		for (r = 0; r < asm_opt.number_of_round; ++r) {
 			clear_opt(&asm_opt, r); // this update asm_opt.roundID and a few other fields
 			ha_overlap_and_correct(r);
-			fprintf(stderr, "[M::%s::%.3f*%.2f] ==> corrected reads for round %d\n", __func__, yak_realtime(), yak_cputime() / yak_realtime(), r + 1);
+			fprintf(stderr, "[M::%s::%.3f*%.2f@%.3fGB] ==> corrected reads for round %d\n", __func__, yak_realtime(),
+					yak_cpu_usage(), yak_peakrss_in_gb(), r + 1);
 			fprintf(stderr, "[M::%s] # bases: %lld; # corrected bases: %lld; # recorrected bases: %lld\n", __func__,
 					asm_opt.num_bases, asm_opt.num_corrected_bases, asm_opt.num_recorrected_bases);
+			fprintf(stderr, "[M::%s] size of buffer: %.3fGB\n", __func__, asm_opt.mem_buf / 1073741824.0);
 		}
 		Output_corrected_reads();
-		fprintf(stderr, "[M::%s::%.3f*%.2f] ==> written corrected reads to disk\n", __func__, yak_realtime(), yak_cputime() / yak_realtime());
+		fprintf(stderr, "[M::%s::%.3f*%.2f] ==> written corrected reads to disk\n", __func__, yak_realtime(), yak_cpu_usage());
 		// overlap between corrected reads
 		clear_opt(&asm_opt, asm_opt.number_of_round);
 		ha_overlap_final();
-		fprintf(stderr, "[M::%s::%.3f*%.2f] ==> found overlaps for the final round\n", __func__, yak_realtime(), yak_cputime() / yak_realtime());
+		fprintf(stderr, "[M::%s::%.3f*%.2f@%.3fGB] ==> found overlaps for the final round\n", __func__, yak_realtime(),
+				yak_cpu_usage(), yak_peakrss_in_gb());
 		ha_print_ovlp_stat(R_INF.paf, R_INF.reverse_paf, R_INF.total_reads);
 		ha_ft_destroy(ha_flt_tab);
 		Output_PAF();

CommandLines.cpp

@@ -105,6 +105,7 @@ void clear_opt(hifiasm_opt_t* asm_opt, int round)
     asm_opt->num_bases = 0;
     asm_opt->num_corrected_bases = 0;
     asm_opt->num_recorrected_bases = 0;
+	asm_opt->mem_buf = 0;
     asm_opt->roundID = round;
 }
 

CommandLines.h

@@ -3,7 +3,7 @@
 
 #include <pthread.h>
 
-#define HA_VERSION "0.3.0-dirty-r191"
+#define HA_VERSION "0.3.0-dirty-r192"
 
 #define VERBOSE 0
 #define VERBOSE_GFA 0
@@ -48,6 +48,7 @@ typedef struct {
     long long num_bases;
     long long num_corrected_bases;
     long long num_recorrected_bases;
+	long long mem_buf;
     long long coverage;
 } hifiasm_opt_t;
 

Hash_Table.cpp

@@ -389,9 +389,9 @@ int32_t ha_chain_lis_core(k_mer_hit *a, int32_t n_a, Chain_Data *dp, int32_t min
 		gap_rate = (double)tot_indel / tot_len;
 		score -= (int)(gap_rate * score * bw_thres);
 		dp->score[i] = dp->score[i-1] + score;
-        dp->pre[i] = i - 1;
-        dp->indels[i] = tot_indel;
-        dp->self_length[i] = tot_len;
+		dp->pre[i] = i - 1;
+		dp->indels[i] = tot_indel;
+		dp->self_length[i] = tot_len;
 	}
 	if (i < m) return -1;
 	for (i = 0; i < m; ++i) a[i] = a[b[i]];

anchor.cpp

@@ -43,6 +43,11 @@ void ha_abuf_destroy(ha_abuf_t *ab)
 	free(ab->seed); free(ab->a); free(ab->mz.a); free(ab);
 }
 
+uint64_t ha_abuf_mem(const ha_abuf_t *ab)
+{
+	return ab->m_a * sizeof(anchor1_t) + ab->mz.m * (sizeof(ha_mz1_t) + sizeof(seed1_t)) + sizeof(ha_abuf_t);
+}
+
 void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_region_alloc *overlap_list, Candidates_list *cl, double bw_thres, int max_n_chain, int keep_whole_chain)
 {
 	extern void *ha_flt_tab;

htab.cpp

@@ -657,7 +657,7 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip
 		free(s->buf);
 		#if 0
 		fprintf(stderr, "[M::%s::%.3f*%.2f] processed %ld sequences; %ld %s in the hash table\n", __func__,
-				yak_realtime(), yak_cputime() / yak_realtime(), (long)s->n_seq0 + s->n_seq,
+				yak_realtime(), yak_cpu_usage(), (long)s->n_seq0 + s->n_seq,
 				(long)(p->pt? p->pt->tot_pos : p->ct->tot), p->pt? "positions" : "distinct k-mers");
 		#endif
 		free(s);
@@ -786,8 +786,8 @@ void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs)
 	ha_ct_shrink(h, cutoff, YAK_MAX_COUNT, asm_opt->thread_num);
 	flt_tab = gen_hh(h);
 	ha_ct_destroy(h);
-	fprintf(stderr, "[M::%s::%.3f*%.2f] ==> filtered out %ld k-mers occurring %d or more times\n", __func__,
-			yak_realtime(), yak_cputime() / yak_realtime(), (long)kh_size(flt_tab), cutoff);
+	fprintf(stderr, "[M::%s::%.3f*%.2f@%.3fGB] ==> filtered out %ld k-mers occurring %d or more times\n", __func__,
+			yak_realtime(), yak_cpu_usage(), yak_peakrss_in_gb(), (long)kh_size(flt_tab), cutoff);
 	return (void*)flt_tab;
 }
 
@@ -808,7 +808,7 @@ ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_f
 	}
 	ct = ha_count(asm_opt, HAF_COUNT_EXACT|extra_flag1, NULL, flt_tab, rs);
 	fprintf(stderr, "[M::%s::%.3f*%.2f] ==> counted %ld distinct minimizer k-mers\n", __func__,
-			yak_realtime(), yak_cputime() / yak_realtime(), (long)ct->tot);
+			yak_realtime(), yak_cpu_usage(), (long)ct->tot);
 	ha_ct_hist(ct, cnt, asm_opt->thread_num);
 	fprintf(stderr, "[M::%s] count[%d] = %ld (for sanity check)\n", __func__, YAK_MAX_COUNT, (long)cnt[YAK_MAX_COUNT]);
 	peak_hom = ha_analyze_count(YAK_N_COUNTS, cnt, &peak_het);
@@ -827,6 +827,6 @@ ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_f
 	assert((uint64_t)tot_cnt == pt->tot_pos);
 	//ha_pt_sort(pt, asm_opt->thread_num);
 	fprintf(stderr, "[M::%s::%.3f*%.2f] ==> indexed %ld positions\n", __func__,
-			yak_realtime(), yak_cputime() / yak_realtime(), (long)pt->tot_pos);
+			yak_realtime(), yak_cpu_usage(), (long)pt->tot_pos);
 	return pt;
 }

htab.h

@@ -36,11 +36,14 @@ const ha_idxpos_t *ha_pt_get(const ha_pt_t *h, uint64_t hash, int *n);
 
 ha_abuf_t *ha_abuf_init(void);
 void ha_abuf_destroy(ha_abuf_t *ab);
+uint64_t ha_abuf_mem(const ha_abuf_t *ab);
 
 double yak_cputime(void);
 void yak_reset_realtime(void);
 double yak_realtime(void);
 long yak_peakrss(void);
+double yak_peakrss_in_gb(void);
+double yak_cpu_usage(void);
 
 void trio_partition(void);
 

main.cpp

@@ -18,6 +18,6 @@ int main(int argc, char *argv[])
 	fprintf(stderr, "[M::%s] CMD:", __func__);
 	for (i = 0; i < argc; ++i)
 		fprintf(stderr, " %s", argv[i]);
-	fprintf(stderr, "\n[M::%s] Real time: %.3f sec; CPU: %.3f sec; Peak RSS: %.3f GB\n", __func__, yak_realtime(), yak_cputime(), yak_peakrss() / 1024.0 / 1024.0 / 1024.0);
+	fprintf(stderr, "\n[M::%s] Real time: %.3f sec; CPU: %.3f sec; Peak RSS: %.3f GB\n", __func__, yak_realtime(), yak_cputime(), yak_peakrss_in_gb());
     return ret;
 }

sys.cpp

@@ -41,3 +41,13 @@ long yak_peakrss(void)
 	return r.ru_maxrss;
 #endif
 }
+
+double yak_peakrss_in_gb(void)
+{
+	return yak_peakrss() / 1073741824.0;
+}
+
+double yak_cpu_usage(void)
+{
+	return (yak_cputime() + 1e-9) / (yak_realtime() + 1e-9);
+}