Context Navigation

source: branches/profile/ARBDB/adseqcompr.cxx

Visit:

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

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

Download in other formats: