Context Navigation

source: tags/ms_r16q2/ARBDB/adseqcompr.cxx

Visit:

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

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

Download in other formats: