Context Navigation

source: tags/arb-6.0-rc1/ARBDB/adseqcompr.cxx

Visit:

Last change on this file was 11838, checked in by westram, 11 years ago
code visibility made functions static moved functions into UNIT_TESTS section
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 39.7 KB

Line
1	// =============================================================== //
2	// //
3	// File : adseqcompr.cxx //
4	// Purpose : //
5	// //
6	// Institute of Microbiology (Technical University Munich) //
7	// http://www.arb-home.de/ //
8	// //
9	// =============================================================== //
10
11	#include <arbdbt.h>
12	#include <arb_progress.h>
13	#include <arb_file.h>
14	#include <arb_misc.h>
15	#include <arb_diff.h>
16
17	#include "gb_key.h"
18	#include <climits>
19
20	// --------------------------------------------------------------------------------
21
22	#define MAX_SEQUENCE_PER_MASTER 50 // was 18 till May 2008
23
24	#if defined(DEBUG)
25	// don't do optimize, only create tree and save to DB
26	// #define SAVE_COMPRESSION_TREE_TO_DB
27	#endif // DEBUG
28
29	// --------------------------------------------------------------------------------
30
31	struct CompressionTree : public GBT_TREE {
32	// members initialized by init_indices_and_count_sons
33	int index; // master(inner nodes) or sequence(leaf nodes) index
34	int sons; // sons with sequence or masters (in subtree)
35
36	~CompressionTree() OVERRIDE {}
37
38	DEFINE_SIMPLE_TREE_RELATIVES_ACCESSORS(CompressionTree);
39	};
40
41	struct Consensus {
42	int len;
43	char used[256];
44	unsigned char *con[256];
45	};
46
47	struct Sequence {
48	GBENTRY *gb_seq;
49	int master;
50	};
51
52	struct MasterSequence {
53	GBENTRY *gb_mas;
54	int master;
55	};
56
57	// --------------------------------------------------------------------------------
58
59	static Consensus *g_b_new_Consensus(long len) {
60	Consensus gcon = (Consensus )GB_calloc(sizeof(*gcon), 1);
61	unsigned char data = (unsigned char )GB_calloc(sizeof(char)*256, len);
62
63	gcon->len = len;
64
65	for (int i=0; i<256; i++) {
66	gcon->con[i] = data + len*i;
67	}
68	return gcon;
69	}
70
71
72	static void g_b_delete_Consensus(Consensus *gcon) {
73	free(gcon->con[0]);
74	free(gcon);
75	}
76
77
78	static void g_b_Consensus_add(Consensus gcon, unsigned char seq, long seq_len) {
79	const int max_priority = 255/MAX_SEQUENCE_PER_MASTER; // No overflow possible
80	gb_assert(max_priority >= 1);
81
82	if (seq_len > gcon->len) seq_len = gcon->len;
83
84	// Search for runs
85	unsigned char *s = seq;
86	int last = 0;
87	int i;
88	int li;
89	int c;
90
91	for (li = i = 0; i < seq_len; i++) {
92	c = *(s++);
93	if (c == last) {
94	continue;
95	}
96	else {
97	inc_hits :
98	int eq_count = i-li;
99	gcon->used[c] = 1;
100	unsigned char *p = gcon->con[last];
101	last = c;
102
103	if (eq_count <= GB_RUNLENGTH_SIZE) {
104	c = max_priority;
105	while (li < i) p[li++] += c;
106	}
107	else {
108	c = max_priority * (GB_RUNLENGTH_SIZE) / eq_count;
109	if (c) {
110	while (li < i) p[li++] += c;
111	}
112	else {
113	while (li < i) p[li++] \|= 1;
114	}
115	}
116	}
117	}
118	if (li<seq_len) {
119	c = last;
120	i = seq_len;
121	goto inc_hits;
122	}
123	}
124
125	static char g_b_Consensus_get_sequence(Consensus gcon) {
126	int pos;
127	unsigned char *s;
128	unsigned char max = (unsigned char )GB_calloc(sizeof(char), gcon->len);
129	int c;
130	char seq = (char )GB_calloc(sizeof(char), gcon->len+1);
131
132	memset(seq, '@', gcon->len);
133
134	for (c = 1; c<256; c++) { // Find maximum frequency of non run
135	if (!gcon->used[c]) continue;
136	s = gcon->con[c];
137	for (pos = 0; pos<gcon->len; pos++) {
138	if (*s > max[pos]) {
139	max[pos] = *s;
140	seq[pos] = c;
141	}
142	s++;
143	}
144	}
145	free(max);
146	return seq;
147	}
148
149
150	static int g_b_count_leafs(CompressionTree *node) {
151	if (node->is_leaf) return 1;
152	node->gb_node = 0;
153	return (g_b_count_leafs(node->get_leftson()) + g_b_count_leafs(node->get_rightson()));
154	}
155
156	static void g_b_put_sequences_in_container(CompressionTree ctree, Sequence seqs, MasterSequence *masters, Consensus gcon) {
157	if (ctree->is_leaf) {
158	if (ctree->index >= 0) {
159	GB_CSTR data = GB_read_char_pntr(seqs[ctree->index].gb_seq);
160	long len = GB_read_string_count(seqs[ctree->index].gb_seq);
161	g_b_Consensus_add(gcon, (unsigned char *)data, len);
162	}
163	}
164	else if (ctree->index<0) {
165	g_b_put_sequences_in_container(ctree->get_leftson(), seqs, masters, gcon);
166	g_b_put_sequences_in_container(ctree->get_rightson(), seqs, masters, gcon);
167	}
168	else {
169	GB_CSTR data = GB_read_char_pntr(masters[ctree->index]->gb_mas);
170	long len = GB_read_string_count(masters[ctree->index]->gb_mas);
171	g_b_Consensus_add(gcon, (unsigned char *)data, len);
172	}
173	}
174
175	static void g_b_create_master(CompressionTree node, Sequence seqs, MasterSequence *masters, int my_master, const char ali_name, long seq_len, arb_progress& progress) {
176	if (node->is_leaf) {
177	if (node->index >= 0) {
178	GBDATA *gb_data = GBT_read_sequence(node->gb_node, ali_name);
179
180	seqs[node->index].gb_seq = gb_data->as_entry();
181	seqs[node->index].master = my_master;
182	}
183	}
184	else {
185	if (progress.aborted()) return;
186
187	if (node->index>=0) {
188	masters[node->index]->master = my_master;
189	my_master = node->index;
190	}
191	g_b_create_master(node->get_leftson(), seqs, masters, my_master, ali_name, seq_len, progress);
192	g_b_create_master(node->get_rightson(), seqs, masters, my_master, ali_name, seq_len, progress);
193	if (node->index>=0 && !progress.aborted()) { // build me
194	char *data;
195	Consensus *gcon = g_b_new_Consensus(seq_len);
196
197	g_b_put_sequences_in_container(node->get_leftson(), seqs, masters, gcon);
198	g_b_put_sequences_in_container(node->get_rightson(), seqs, masters, gcon);
199
200	data = g_b_Consensus_get_sequence(gcon);
201
202	GB_write_string(masters[node->index]->gb_mas, data);
203	GB_write_security_write(masters[node->index]->gb_mas, 7);
204
205	g_b_delete_Consensus(gcon);
206	free(data);
207
208	++progress;
209	}
210	}
211	}
212
213	// -------------------------------------
214	// distribute master sequences
215
216	static void subtract_sons_from_tree(CompressionTree *node, int subtract) {
217	while (node) {
218	node->sons -= subtract;
219	node = node->get_father();
220	}
221	}
222
223	static int set_masters_with_sons(CompressionTree node, int wantedSons, int mcount) {
224	if (!node->is_leaf) {
225	if (node->sons == wantedSons) {
226	// insert new master
227	gb_assert(node->index == -1);
228	node->index = *mcount;
229	(*mcount)++;
230
231	subtract_sons_from_tree(node->get_father(), node->sons-1);
232	node->sons = 1;
233	}
234	else if (node->sons>wantedSons) {
235	int lMax = set_masters_with_sons(node->get_leftson(), wantedSons, mcount);
236	int rMax = set_masters_with_sons(node->get_rightson(), wantedSons, mcount);
237
238	int maxSons = lMax<rMax ? rMax : lMax;
239	if (node->sons <= MAX_SEQUENCE_PER_MASTER && node->sons>maxSons) {
240	maxSons = node->sons;
241	}
242	return maxSons;
243	}
244	}
245	return node->sons <= MAX_SEQUENCE_PER_MASTER ? node->sons : 0;
246	}
247
248	static int maxCompressionSteps(CompressionTree *node) {
249	if (node->is_leaf) {
250	return 0;
251	}
252
253	int left = maxCompressionSteps(node->get_leftson());
254	int right = maxCompressionSteps(node->get_rightson());
255
256	#if defined(SAVE_COMPRESSION_TREE_TO_DB)
257	freenull(node->name);
258	if (node->index2 != -1) {
259	node->name = GBS_global_string_copy("master_%03i", node->index2);
260	}
261	#endif // SAVE_COMPRESSION_TREE_TO_DB
262
263	return (left>right ? left : right) + (node->index == -1 ? 0 : 1);
264	}
265
266	static int init_indices_and_count_sons(CompressionTree node, int scount, const char *ali_name) {
267	if (node->is_leaf) {
268	if (node->gb_node == 0 \|\| !GBT_read_sequence(node->gb_node, (char *)ali_name)) {
269	node->index = -1;
270	node->sons = 0;
271	}
272	else {
273	node->index = *scount;
274	node->sons = 1;
275	(*scount)++;
276	}
277	}
278	else {
279	node->index = -1;
280	node->sons =
281	init_indices_and_count_sons(node->get_leftson(), scount, ali_name) +
282	init_indices_and_count_sons(node->get_rightson(), scount, ali_name);
283	}
284	return node->sons;
285	}
286
287	static void distribute_masters(CompressionTree tree, int mcount, int *max_masters) {
288	int wantedSons = MAX_SEQUENCE_PER_MASTER;
289	while (wantedSons >= 2) {
290	int maxSons = set_masters_with_sons(tree, wantedSons, mcount);
291	wantedSons = maxSons;
292	}
293	gb_assert(tree->sons == 1);
294
295	gb_assert(tree->index != -1);
296	*max_masters = maxCompressionSteps(tree);
297	}
298
299	// --------------------------------------------------------------------------------
300
301	#define MAX_NUMBER 0x7fffffff
302	STATIC_ASSERT(MAX_NUMBER <= INT_MAX); // ensure 32-bit-version compatibility!
303
304	inline int g_b_read_number2(const unsigned char*& s) {
305	int result;
306	unsigned c0 = *s++;
307	if (c0 & 0x80) {
308	unsigned c1 = *s++;
309	if (c0 & 0x40) {
310	unsigned c2 = *s++;
311	if (c0 & 0x20) {
312	unsigned c3 = *s++;
313	if (c0 & 0x10) {
314	unsigned c4 = *s++;
315	result = c4 \| (c3<<8) \| (c2<<16) \| (c1<<24);
316	}
317	else {
318	result = c3 \| (c2<<8) \| (c1<<16) \| ((c0 & 0x0f)<<24);
319	}
320	}
321	else {
322	result = c2 \| (c1<<8) \| ((c0 & 0x1f)<<16);
323	}
324	}
325	else {
326	result = c1 \| ((c0 & 0x3f)<<8);
327	}
328	}
329	else {
330	result = c0;
331	}
332	gb_assert(result >= 0 && result <= MAX_NUMBER);
333	return result;
334	}
335
336	inline void g_b_put_number2(int i, unsigned char*& s) {
337	gb_assert(i >= 0 && i <= MAX_NUMBER);
338
339	if (i< 0x80) {
340	*s++ = i;
341	}
342	else {
343	int j;
344	if (i<0x4000) {
345	j = (i>>8) \| 0x80; *s++ = j;
346	*s++ = i;
347	}
348	else if (i<0x200000) {
349	j = (i>>16) \| 0xC0; *s++ = j;
350	j = (i>>8); *s++ = j;
351	*s++ = i;
352	}
353	else if (i<0x10000000) {
354	j = (i>>24) \| 0xE0; *s++ = j;
355	j = (i>>16); *s++ = j;
356	j = (i>>8); *s++ = j;
357	*s++ = i;
358	}
359	else {
360	*s++ = 0xF0;
361	j = (i>>24); *s++ = j;
362	j = (i>>16); *s++ = j;
363	j = (i>>8); *s++ = j;
364	*s++ = i;
365	}
366	}
367	}
368
369	// --------------------------------------------------------------------------------
370
371	#ifdef UNIT_TESTS
372	#ifndef TEST_UNIT_H
373	#include <test_unit.h>
374	#endif
375
376	static arb_test::match_expectation put_read_num_using_bytes(int num_written, int bytes_expected, unsigned char *buffer_expected = NULL) {
377	const int BUFSIZE = 6;
378	unsigned char buffer[BUFSIZE];
379
380	using namespace arb_test;
381
382	unsigned char INIT = 0xaa;
383	memset(buffer, INIT, BUFSIZE);
384
385	expectation_group expected;
386
387	{
388	unsigned char *bp = buffer;
389	g_b_put_number2(num_written, bp);
390
391	size_t bytes_written = bp-buffer;
392	expected.add(that(bytes_written).is_equal_to(bytes_expected));
393
394	if (buffer_expected) {
395	expected.add(that(arb_test::memory_is_equal(buffer, buffer_expected, bytes_expected)).is_equal_to(true));
396	}
397	}
398	{
399	const unsigned char *bp = buffer;
400
401	int num_read = g_b_read_number2(bp);
402	expected.add(that(num_read).is_equal_to(num_written));
403
404	size_t bytes_read = bp-buffer;
405	expected.add(that(bytes_read).is_equal_to(bytes_expected));
406	}
407
408	expected.add(that(buffer[bytes_expected]).is_equal_to(INIT));
409
410	return all().ofgroup(expected);
411	}
412
413	#define TEST_PUT_READ_NUMBER(num,expect_bytes) TEST_EXPECTATION(put_read_num_using_bytes(num, expect_bytes))
414	#define TEST_PUT_READ_NUMBER__BROKEN(num,expect_bytes) TEST_EXPECTATION__BROKEN(put_read_num_using_bytes(num, expect_bytes))
415
416	#define TEST_PUT_NUMBER_BINARY1(num, byte1) do { \
417	unsigned char buf[1]; \
418	buf[0] = byte1; \
419	TEST_EXPECTATION(put_read_num_using_bytes(num, 1, buf)); \
420	} while(0)
421
422	#define TEST_PUT_NUMBER_BINARY2(num, byte1, byte2) do { \
423	unsigned char buf[2]; \
424	buf[0] = byte1; \
425	buf[1] = byte2; \
426	TEST_EXPECTATION(put_read_num_using_bytes(num, 2, buf)); \
427	} while(0)
428
429	#define TEST_PUT_NUMBER_BINARY3(num, byte1, byte2, byte3) do { \
430	unsigned char buf[3]; \
431	buf[0] = byte1; \
432	buf[1] = byte2; \
433	buf[2] = byte3; \
434	TEST_EXPECTATION(put_read_num_using_bytes(num, 3, buf)); \
435	} while(0)
436
437	#define TEST_PUT_NUMBER_BINARY4(num, byte1, byte2, byte3, byte4) do { \
438	unsigned char buf[4]; \
439	buf[0] = byte1; \
440	buf[1] = byte2; \
441	buf[2] = byte3; \
442	buf[3] = byte4; \
443	TEST_EXPECTATION(put_read_num_using_bytes(num, 4, buf)); \
444	} while(0)
445
446	#define TEST_PUT_NUMBER_BINARY5(num, byte1, byte2, byte3, byte4, byte5) do { \
447	unsigned char buf[5]; \
448	buf[0] = byte1; \
449	buf[1] = byte2; \
450	buf[2] = byte3; \
451	buf[3] = byte4; \
452	buf[4] = byte5; \
453	TEST_EXPECTATION(put_read_num_using_bytes(num, 5, buf)); \
454	} while(0)
455
456	void TEST_put_read_number() {
457	// test that put and read are compatible:
458	TEST_PUT_READ_NUMBER(0x0, 1);
459
460	TEST_PUT_READ_NUMBER(0x7f, 1);
461	TEST_PUT_READ_NUMBER(0x80, 2);
462
463	TEST_PUT_READ_NUMBER(0x3fff, 2);
464	TEST_PUT_READ_NUMBER(0x4000, 3);
465
466	TEST_PUT_READ_NUMBER(0x1fffff, 3);
467	TEST_PUT_READ_NUMBER(0x200000, 4);
468
469	TEST_PUT_READ_NUMBER(0xfffffff, 4);
470
471	TEST_PUT_READ_NUMBER(0x10000000, 5);
472	TEST_PUT_READ_NUMBER(0x7fffffff, 5);
473
474	// test binary compatibility:
475	// (code affects DB content, cannot be changed)
476
477	TEST_PUT_NUMBER_BINARY1(0x0, 0x00);
478	TEST_PUT_NUMBER_BINARY1(0x7f, 0x7f);
479
480	TEST_PUT_NUMBER_BINARY2(0x80, 0x80, 0x80);
481	TEST_PUT_NUMBER_BINARY2(0x81, 0x80, 0x81);
482	TEST_PUT_NUMBER_BINARY2(0x3fff, 0xbf, 0xff);
483
484	TEST_PUT_NUMBER_BINARY3(0x4000, 0xc0, 0x40, 0x00);
485	TEST_PUT_NUMBER_BINARY3(0x1fffff, 0xdf, 0xff, 0xff);
486
487	TEST_PUT_NUMBER_BINARY4(0x200000, 0xe0, 0x20, 0x00, 0x00);
488	TEST_PUT_NUMBER_BINARY4(0xfffffff, 0xef, 0xff, 0xff, 0xff);
489
490	TEST_PUT_NUMBER_BINARY5(0x10000000, 0xf0, 0x10, 0x00, 0x00, 0x00);
491	TEST_PUT_NUMBER_BINARY5(0x7fffffff, 0xf0, 0x7f, 0xff, 0xff, 0xff);
492	}
493
494	#endif // UNIT_TESTS
495
496	// --------------------------------------------------------------------------------
497
498
499	static char gb_compress_seq_by_master(const char master, size_t master_len, int master_index,
500	GBQUARK q, const char *seq, size_t seq_len,
501	size_t *memsize, int old_flag) {
502	unsigned char *buffer;
503	int rest = 0;
504	unsigned char *d;
505	int i, cs, cm;
506	int last;
507	long len = seq_len;
508
509	d = buffer = (unsigned char *)GB_give_other_buffer(seq, seq_len);
510
511	if (seq_len > master_len) {
512	rest = seq_len - master_len;
513	len = master_len;
514	}
515
516	last = -1000; // Convert Sequence relative to Master
517	for (i = len; i>0; i--) {
518	cm = *(master++);
519	cs = *(seq++);
520	if (cm==cs && cs != last) {
521	*(d++) = 0;
522	last = 1000;
523	}
524	else {
525	*(d++) = cs;
526	last = cs;
527	}
528	}
529	for (i = rest; i>0; i--) {
530	(d++) = (seq++);
531	}
532
533	{ // Append run length compression method
534	unsigned char *buffer2;
535	unsigned char *dest2;
536	buffer2 = dest2 = (unsigned char )GB_give_other_buffer((char )buffer, seq_len+100);
537	*(dest2++) = GB_COMPRESSION_SEQUENCE \| old_flag;
538
539	g_b_put_number2(master_index, dest2); // Tags
540	g_b_put_number2(q, dest2);
541
542	gb_compress_equal_bytes_2((char )buffer, seq_len, memsize, (char )dest2); // append runlength compressed sequences to tags
543
544	memsize = memsize + (dest2-buffer2);
545	return (char *)buffer2;
546	}
547	}
548
549	static char gb_compress_sequence_by_master(GBDATA gbd, const char *master, size_t master_len, int master_index,
550	GBQUARK q, const char seq, size_t seq_len, size_t memsize)
551	{
552	size_t size;
553	char *is = gb_compress_seq_by_master(master, master_len, master_index, q, seq, seq_len, &size, GB_COMPRESSION_LAST);
554	char *res = gb_compress_data(gbd, 0, is, size, memsize, ~(GB_COMPRESSION_DICTIONARY\|GB_COMPRESSION_SORTBYTES\|GB_COMPRESSION_RUNLENGTH), true);
555	return res;
556	}
557
558	static GB_ERROR compress_sequence_tree(GBCONTAINER gb_main, CompressionTree tree, const char *ali_name) {
559	GB_ERROR error = 0;
560	long ali_len = GBT_get_alignment_len(gb_main, ali_name);
561	int main_clock = GB_read_clock(gb_main);
562
563	GB_ERROR warning = NULL;
564
565	if (ali_len<0) {
566	warning = GBS_global_string("Skipping alignment '%s' (not a valid alignment; len=%li).", ali_name, ali_len);
567	}
568	else {
569	int leafcount = g_b_count_leafs(tree);
570	if (!leafcount) {
571	error = "Tree is empty";
572	}
573	else {
574	arb_progress tree_progress("Compressing sequences", 4);
575
576	// Distribute masters in tree
577	int mastercount = 0;
578	int max_compSteps = 0; // in one branch
579	int seqcount = 0;
580
581	init_indices_and_count_sons(tree, &seqcount, ali_name);
582	if (!seqcount) {
583	warning = GBS_global_string("Tree contains no sequences with data in '%s'\n"
584	"Skipping compression for this alignment",
585	ali_name);
586	}
587	else {
588	distribute_masters(tree, &mastercount, &max_compSteps);
589
590	#if defined(SAVE_COMPRESSION_TREE_TO_DB)
591	{
592	error = GBT_link_tree(tree, gb_main, 0, NULL, NULL);
593	if (!error) error = GBT_write_tree(gb_main, 0, "tree_compression_new", tree);
594	GB_information("Only generated compression tree (do NOT save DB anymore)");
595	return error;
596	}
597	#endif // SAVE_COMPRESSION_TREE_TO_DB
598
599	// detect degenerated trees
600	{
601	int min_masters = ((seqcount-1)/MAX_SEQUENCE_PER_MASTER)+1;
602	int min_compSteps = 1;
603	{
604	int m = min_masters;
605	while (m>1) {
606	m = ((m-1)/MAX_SEQUENCE_PER_MASTER)+1;
607	min_masters += m;
608	min_compSteps++;
609	}
610	}
611
612	int acceptable_masters = (3*min_masters)/2; // accept 50% overhead
613	int acceptable_compSteps = 11*min_compSteps; // accept 1000% overhead
614
615	if (mastercount>acceptable_masters \|\| max_compSteps>acceptable_compSteps) {
616	GB_warningf("Tree is ill-suited for compression (cause of deep branches)\n"
617	" Used tree Optimal tree Overhead\n"
618	"Compression steps %5i %5i %4i%% (speed)\n"
619	"Master sequences %5i %5i %4i%% (size)\n"
620	"If you like to restart with a better tree,\n"
621	"press 'Abort' to stop compression",
622	max_compSteps, min_compSteps, (100*max_compSteps)/min_compSteps-100,
623	mastercount, min_masters, (100*mastercount)/min_masters-100);
624	}
625	}
626
627	gb_assert(mastercount>0);
628	}
629
630	if (!warning) {
631	GBCONTAINER *gb_master_ali = 0;
632	GBDATA *old_gb_master_ali = 0;
633	Sequence *seqs = 0;
634	GB_MAIN_TYPE *Main = GB_MAIN(gb_main);
635	GBQUARK ali_quark = gb_find_or_create_quark(Main, ali_name);
636	unsigned long long sumorg = 0;
637	unsigned long long sumold = 0;
638	unsigned long long sumnew = 0;
639	MasterSequence masters = (MasterSequence )GB_calloc(sizeof(*masters), leafcount);
640	int si;
641
642	{
643	char *masterfoldername = GBS_global_string_copy("%s/@master_data/@%s", GB_SYSTEM_FOLDER, ali_name);
644	old_gb_master_ali = GB_search(gb_main, masterfoldername, GB_FIND)->as_container();
645	free(masterfoldername);
646	}
647
648	// create masters
649	if (!error) {
650	{
651	char *master_data_name = GBS_global_string_copy("%s/@master_data", GB_SYSTEM_FOLDER);
652	char *master_name = GBS_global_string_copy("@%s", ali_name);
653
654	GBCONTAINER *gb_master_data = gb_search(gb_main, master_data_name, GB_CREATE_CONTAINER, 1)->as_container();
655
656	// create a master container, the old is deleted as soon as all sequences are compressed by the new method
657	gb_master_ali = gb_create_container(gb_master_data, master_name);
658	GB_write_security_delete(gb_master_ali, 7);
659
660	free(master_name);
661	free(master_data_name);
662	}
663	for (si = 0; si<mastercount; si++) {
664	masters[si] = (MasterSequence *)GB_calloc(sizeof(MasterSequence), 1);
665	masters[si]->gb_mas = gb_create(gb_master_ali, "@master", GB_STRING);
666	}
667	seqs = (Sequence )GB_calloc(sizeof(seqs), leafcount);
668
669	if (!error) {
670	arb_progress progress("Building master sequences", mastercount);
671	g_b_create_master(tree, seqs, masters, -1, ali_name, ali_len, progress);
672
673	error = progress.error_if_aborted();
674	}
675	}
676	tree_progress.inc_and_check_user_abort(error);
677
678	// Compress sequences in tree
679	if (!error) {
680	arb_progress progress("Compressing sequences in tree", seqcount);
681
682	for (si=0; si<seqcount && !error; si++) {
683	int mi = seqs[si].master;
684	MasterSequence *master = masters[mi];
685	GBDATA *gbd = seqs[si].gb_seq;
686
687	if (GB_read_clock(gbd) >= main_clock) {
688	GB_warning("A species seems to be more than once in the tree");
689	}
690	else {
691	char *seq = GB_read_string(gbd);
692	int seq_len = GB_read_string_count(gbd);
693	long sizen = GB_read_memuse(gbd);
694	char *seqm = GB_read_string(master->gb_mas);
695	int master_len = GB_read_string_count(master->gb_mas);
696	size_t sizes;
697	char *ss = gb_compress_sequence_by_master(gbd, seqm, master_len, mi, ali_quark, seq, seq_len, &sizes);
698
699	gb_write_compressed_pntr(gbd->as_entry(), ss, sizes, seq_len);
700	sizes = GB_read_memuse(gbd); // check real usage
701
702	sumnew += sizes;
703	sumold += sizen;
704	sumorg += seq_len;
705
706	free(seqm);
707	free(seq);
708	}
709
710	progress.inc_and_check_user_abort(error);
711	}
712	}
713	tree_progress.inc_and_check_user_abort(error);
714
715	// Compress rest of sequences
716	if (!error) {
717	int pass; // pass 1 : count species to compress, pass 2 : compress species
718	int speciesNotInTree = 0;
719
720	SmartPtr<arb_progress> progress;
721
722	for (pass = 1; pass <= 2; ++pass) {
723	GBDATA *gb_species_data = GBT_get_species_data(gb_main);
724	GBDATA *gb_species;
725	int count = 0;
726
727	for (gb_species = GBT_first_species_rel_species_data(gb_species_data);
728	gb_species;
729	gb_species = GBT_next_species(gb_species))
730	{
731	GBDATA *gbd = GBT_read_sequence(gb_species, ali_name);
732
733	if (!gbd) continue;
734	if (GB_read_clock(gbd) >= main_clock) continue; // Compress only those which are not compressed by masters
735	count++;
736	if (pass == 2) {
737	char *data = GB_read_string(gbd);
738	long seq_len = GB_read_string_count(gbd);
739	long size = GB_read_memuse(gbd);
740
741	GB_write_string(gbd, "");
742	GB_write_string(gbd, data);
743	free(data);
744
745	sumold += size;
746
747	size = GB_read_memuse(gbd);
748	sumnew += size;
749	sumorg += seq_len;
750
751	progress->inc_and_check_user_abort(error);
752	}
753	}
754	if (pass == 1) {
755	speciesNotInTree = count;
756	if (speciesNotInTree>0) {
757	progress = new arb_progress("Compressing sequences NOT in tree", speciesNotInTree);
758	}
759	}
760	}
761	}
762	tree_progress.inc_and_check_user_abort(error);
763
764	if (!error) {
765	arb_progress progress("Compressing master-sequences", mastercount);
766
767	// Compress all masters
768	for (si=0; si<mastercount; si++) {
769	int mi = masters[si]->master;
770
771	if (mi>0) { // master available
772	GBDATA *gbd = masters[si]->gb_mas;
773
774	gb_assert(mi>si); // we don't want a recursion, because we cannot uncompress sequence compressed masters, Main->gb_master_data is wrong
775
776	if (gb_read_nr(gbd) != si) { // Check database
777	GB_internal_error("Sequence Compression: Master Index Conflict");
778	error = GB_export_error("Sequence Compression: Master Index Conflict");
779	break;
780	}
781
782	{
783	MasterSequence *master = masters[mi];
784	char *seqm = GB_read_string(master->gb_mas);
785	int master_len = GB_read_string_count(master->gb_mas);
786	char *seq = GB_read_string(gbd);
787	int seq_len = GB_read_string_count(gbd);
788	size_t sizes;
789	char *ss = gb_compress_sequence_by_master(gbd, seqm, master_len, mi, ali_quark, seq, seq_len, &sizes);
790
791	gb_write_compressed_pntr(gbd->as_entry(), ss, sizes, seq_len);
792	sumnew += sizes;
793
794	free(seq);
795	free(seqm);
796	}
797
798	progress.inc_and_check_user_abort(error);
799	}
800	else { // count size of top master
801	GBDATA *gbd = masters[si]->gb_mas;
802	sumnew += GB_read_memuse(gbd);
803
804	progress.inc_and_check_user_abort(error);
805	}
806	}
807
808	// count size of old master data
809	if (!error) {
810	GBDATA *gb_omaster;
811	for (gb_omaster = GB_entry(old_gb_master_ali, "@master");
812	gb_omaster;
813	gb_omaster = GB_nextEntry(gb_omaster))
814	{
815	long size = GB_read_memuse(gb_omaster);
816	sumold += size;
817	}
818	}
819
820	if (!error) {
821	char *sizeOrg = strdup(GBS_readable_size(sumorg, "b"));
822	char *sizeOld = strdup(GBS_readable_size(sumold, "b"));
823	char *sizeNew = strdup(GBS_readable_size(sumnew, "b"));
824
825	GB_warningf("Alignment '%s':\n"
826	" Uncompressed data: %7s\n"
827	" Old compressed data: %7s = %6.2f%%\n"
828	" New compressed data: %7s = %6.2f%%",
829	ali_name, sizeOrg,
830	sizeOld, (100.0*sumold)/sumorg,
831	sizeNew, (100.0*sumnew)/sumorg);
832
833	free(sizeNew);
834	free(sizeOld);
835	free(sizeOrg);
836	}
837	}
838	tree_progress.inc_and_check_user_abort(error);
839
840	if (!error) {
841	if (old_gb_master_ali) error = GB_delete(old_gb_master_ali);
842	Main->keys[ali_quark].gb_master_ali = gb_master_ali;
843	}
844
845	// free data
846	free(seqs);
847	for (si=0; si<mastercount; si++) free(masters[si]);
848	free(masters);
849	}
850	}
851	}
852
853	if (warning) GB_information(warning);
854
855	return error;
856	}
857
858	class CompressionTree_NodeFactory : public TreeNodeFactory {
859	virtual GBT_TREE *makeNode() const OVERRIDE { return new CompressionTree; }
860	};
861
862	GB_ERROR GBT_compress_sequence_tree2(GBDATA gbd, const char tree_name, const char *ali_name) { // goes to header: __ATTR__USERESULT // @@@ rename function
863	// Compress sequences, call only outside a transaction
864	GB_ERROR error = NULL;
865	GB_MAIN_TYPE *Main = GB_MAIN(gbd);
866
867	if (Main->get_transaction_level() > 0) {
868	error = "Compress Sequences called while transaction running";
869	GB_internal_error(error);
870	}
871	else {
872	GBCONTAINER *gb_main = Main->root_container;
873	GB_UNDO_TYPE prev_undo_type = GB_get_requested_undo_type(gb_main);
874
875	error = GB_request_undo_type(gb_main, GB_UNDO_KILL);
876	if (!error) {
877	error = GB_begin_transaction(gb_main);
878	if (!error) {
879	GB_push_my_security(gb_main);
880
881	if (!tree_name \|\| !strlen(tree_name)) {
882	tree_name = GBT_name_of_largest_tree(gb_main);
883	}
884
885	{
886	CompressionTree ctree = DOWNCAST(CompressionTree, GBT_read_tree(gb_main, tree_name, CompressionTree_NodeFactory()));
887	if (!ctree) error = GB_await_error();
888	else {
889	error = GBT_link_tree(ctree, gb_main, false, 0, 0);
890	if (!error) error = compress_sequence_tree(gb_main, ctree, ali_name);
891	delete ctree;
892	}
893	}
894	if (!error) GB_disable_quicksave(gb_main, "Database optimized");
895
896	GB_pop_my_security(gb_main);
897	error = GB_end_transaction(gb_main, error);
898	}
899	ASSERT_NO_ERROR(GB_request_undo_type(gb_main, prev_undo_type));
900	}
901
902	#if defined(SAVE_COMPRESSION_TREE_TO_DB)
903	error = "fake error";
904	#endif // SAVE_COMPRESSION_TREE_TO_DB
905	}
906	return error;
907	}
908
909	#ifdef DEBUG
910
911	void GBT_compression_test(void /dummy_AW_root/, GBDATA gb_main) {
912	GB_ERROR error = GB_begin_transaction(gb_main);
913	char *ali_name = GBT_get_default_alignment(gb_main);
914	char *tree_name = GBT_read_string(gb_main, "focus/tree_name");
915
916	// GBUSE(dummy);
917	if (!ali_name \|\| !tree_name) error = GB_await_error();
918
919	error = GB_end_transaction(gb_main, error);
920
921	if (!error) {
922	printf("Recompression data in alignment '%s' using tree '%s'\n", ali_name, tree_name);
923	error = GBT_compress_sequence_tree2(gb_main, tree_name, ali_name);
924	}
925
926	if (error) GB_warning(error);
927	free(tree_name);
928	free(ali_name);
929	}
930
931	#endif
932
933	// ****************** Decompress Sequences ******************
934
935	static char g_b_uncompress_single_sequence_by_master(const char s, const char master, size_t size, size_t new_size) {
936	const signed char source = (signed char )s;
937	char *dest;
938	const char *m = master;
939	unsigned int c;
940	int j;
941	int i;
942	char *buffer;
943
944	dest = buffer = GB_give_other_buffer((char *)source, size);
945
946	for (i=size; i;) {
947	j = *(source++);
948	if (j>0) { // uncompressed data block
949	if (j>i) j=i;
950	i -= j;
951	for (; j; j--) {
952	c = *(source++);
953	if (!c) c = *m;
954	*(dest++) = c;
955	m++;
956	}
957	}
958	else { // equal bytes compressed
959	if (!j) break; // end symbol
960	if (j == -122) {
961	j = *(source++) & 0xff;
962	j \|= ((*(source++)) << 8) &0xff00;
963	j = -j;
964	}
965	c = *(source++);
966	i += j;
967	if (i<0) {
968	GB_internal_error("Internal Error: Missing end in data");
969	j += -i;
970	i = 0;
971	}
972	if (c==0) {
973	memcpy(dest, m, -j);
974	dest += -j;
975	m += -j;
976	}
977	else {
978	memset(dest, c, -j);
979	dest += -j;
980	m += -j;
981	}
982	}
983	}
984	*(dest++) = 0; // NULL of NULL terminated string
985
986	*new_size = dest-buffer;
987	gb_assert(size == *new_size); // buffer overflow
988
989	return buffer;
990	}
991
992	char gb_uncompress_by_sequence(GBDATA gbd, const char ss, size_t size, GB_ERROR error, size_t *new_size) {
993	char *dest = 0;
994
995	*error = 0;
996	if (!GB_FATHER(gbd)) {
997	*error = "Cannot uncompress this sequence: Sequence has no father";
998	}
999	else {
1000	GB_MAIN_TYPE *Main = GB_MAIN(gbd);
1001	GBDATA *gb_main = Main->gb_main();
1002	char *to_free = GB_check_out_buffer(ss); // Remove 'ss' from memory management, otherwise load_single_key_data() may destroy it
1003	int index;
1004	GBQUARK quark;
1005
1006	{
1007	const unsigned char s = (const unsigned char )ss;
1008
1009	index = g_b_read_number2(s);
1010	quark = g_b_read_number2(s);
1011
1012	ss = (const char *)s;
1013	}
1014
1015	if (!Main->keys[quark].gb_master_ali) {
1016	gb_load_single_key_data(gb_main, quark);
1017	}
1018
1019	if (!Main->keys[quark].gb_master_ali) {
1020	*error = "Cannot uncompress this sequence: Cannot find a master sequence";
1021	}
1022	else {
1023	GBDATA *gb_master = gb_find_by_nr(Main->keys[quark].gb_master_ali, index);
1024	if (gb_master) {
1025	const char *master = GB_read_char_pntr(gb_master); // make sure that this is not a buffer !!!
1026
1027	gb_assert((GB_read_string_count(gb_master)+1) == size); // size mismatch between master and slave
1028	dest = g_b_uncompress_single_sequence_by_master(ss, master, size, new_size);
1029	}
1030	else {
1031	*error = GB_await_error();
1032	}
1033	}
1034	free(to_free);
1035	}
1036
1037	return dest;
1038	}
1039
1040	// --------------------------------------------------------------------------------
1041
1042	#ifdef UNIT_TESTS
1043	#ifndef TEST_UNIT_H
1044	#include <test_unit.h>
1045	#endif
1046
1047	// #define TEST_AUTO_UPDATE // uncomment to auto-update expected result DB
1048
1049	void TEST_SLOW_sequence_compression() {
1050	const char *source = "TEST_nuc.arb";
1051	const char *compressed = "TEST_nuc_seqcompr.arb";
1052	const char *expected = "TEST_nuc_seqcompr_exp.arb";
1053	const char *aliname = "ali_16s";
1054
1055	GB_shell shell;
1056
1057	const int SEQ2COMPARE = 7;
1058	char *seq_exp[SEQ2COMPARE];
1059
1060	{
1061	GBDATA *gb_main;
1062	TEST_EXPECT_RESULT__NOERROREXPORTED(gb_main = GB_open(source, "rw"));
1063
1064	{
1065	GB_transaction ta(gb_main);
1066	int count = 0;
1067
1068	for (GBDATA *gb_species = GBT_first_species(gb_main);
1069	gb_species && count<SEQ2COMPARE;
1070	gb_species = GBT_next_species(gb_species), ++count)
1071	{
1072	GBDATA *gb_seq = GBT_read_sequence(gb_species, aliname);
1073	seq_exp[count] = GB_read_string(gb_seq);
1074	}
1075	}
1076
1077	TEST_EXPECT_NO_ERROR(GBT_compress_sequence_tree2(gb_main, "tree_nuc", aliname));
1078	TEST_EXPECT_NO_ERROR(GB_save_as(gb_main, compressed, "b"));
1079	GB_close(gb_main);
1080	}
1081	#if defined(TEST_AUTO_UPDATE)
1082	TEST_COPY_FILE(compressed, expected);
1083	#endif
1084	TEST_EXPECT_FILES_EQUAL(compressed, expected);
1085
1086	{
1087	GBDATA *gb_main;
1088	TEST_EXPECT_RESULT__NOERROREXPORTED(gb_main = GB_open(compressed, "rw"));
1089	{
1090	GB_transaction ta(gb_main);
1091	int count = 0;
1092
1093	for (GBDATA *gb_species = GBT_first_species(gb_main);
1094	gb_species && count<SEQ2COMPARE;
1095	gb_species = GBT_next_species(gb_species), ++count)
1096	{
1097	GBDATA *gb_seq = GBT_read_sequence(gb_species, aliname);
1098	char *seq = GB_read_string(gb_seq);
1099
1100	TEST_EXPECT_EQUAL(seq, seq_exp[count]);
1101
1102	freenull(seq_exp[count]);
1103	free(seq);
1104	}
1105	}
1106	GB_close(gb_main);
1107	}
1108
1109	TEST_EXPECT_ZERO_OR_SHOW_ERRNO(GB_unlink(compressed));
1110	}
1111
1112	#endif // UNIT_TESTS
1113
1114	// --------------------------------------------------------------------------------
1115

Note: See TracBrowser for help on using the repository browser.

Download in other formats: