Context Navigation

source: branches/stable/SL/FAST_ALIGNER/fast_aligner.cxx

Visit:

Last change on this file was 18665, checked in by westram, 3 years ago
change many WARN_TODO triggered warnings into todo markers.
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 136.4 KB

Line
1	// =============================================================== //
2	// //
3	// File : fast_aligner.cxx //
4	// Purpose : A fast aligner (not a multiple aligner!) //
5	// //
6	// Coded by Ralf Westram (coder@reallysoft.de) in 1998 //
7	// Institute of Microbiology (Technical University Munich) //
8	// http://www.arb-home.de/ //
9	// //
10	// =============================================================== //
11
12
13	#include "fast_aligner.hxx"
14	#include "ClustalV.hxx"
15	#include "seq_search.hxx"
16
17	#include <island_hopping.h>
18
19	#include <awtc_next_neighbours.hxx>
20	#include <awt_sel_boxes.hxx>
21
22	#include <aw_awars.hxx>
23	#include <aw_root.hxx>
24	#include <aw_question.hxx>
25
26	#include <arbdbt.h>
27	#include <ad_cb.h>
28
29	#include <arb_defs.h>
30	#include <arb_progress.h>
31	#include <RangeList.h>
32
33	#include <cctype>
34	#include <cmath>
35	#include <climits>
36
37	#include <list>
38	#include <awt_config_manager.hxx>
39	#include <rootAsWin.h>
40
41	// --------------------------------------------------------------------------------
42
43	#if defined(DEBUG)
44	// #define TRACE_CLUSTAL_DATA
45	// #define TRACE_ISLANDHOPPER_DATA
46	// #define TRACE_COMPRESSED_ALIGNMENT
47	// #define TRACE_RELATIVES
48	#endif // DEBUG
49
50	// --------------------------------------------------------------------------------
51
52	enum FA_report {
53	FA_NO_REPORT, // no report
54	FA_TEMP_REPORT, // report to temporary entries
55	FA_REPORT, // report to resident entries
56	};
57
58	enum FA_range {
59	FA_WHOLE_SEQUENCE, // align whole sequence
60	FA_AROUND_CURSOR, // align xxx positions around current cursor position
61	FA_SELECTED_RANGE, // align selected range
62	FA_SAI_MULTI_RANGE, // align versus multi range define by SAI
63	};
64
65	enum FA_turn {
66	FA_TURN_NEVER, // never try to turn sequence
67	FA_TURN_INTERACTIVE, // try to turn, but query user
68	FA_TURN_ALWAYS, // turn if score is better
69	};
70
71	enum FA_reference {
72	FA_REF_EXPLICIT, // reference sequence explicitly specified
73	FA_REF_CONSENSUS, // use group consensus as reference
74	FA_REF_RELATIVES, // search next relatives by PT server
75	};
76
77	enum FA_alignTarget {
78	FA_CURRENT, // align current species
79	FA_MARKED, // align all marked species
80	FA_SELECTED, // align selected species (= range)
81	};
82
83	enum FA_errorAction {
84	FA_NO_ACTION, // do nothing
85	FA_MARK_FAILED, // mark failed species (unmark rest)
86	FA_MARK_ALIGNED, // mark aligned species (unmark rest)
87	};
88
89	struct AlignParams {
90	FA_report report;
91	bool showGapsMessages; // display messages about missing gaps in master?
92	PosRange range; // range to be aligned
93	};
94
95	struct SearchRelativeParams : virtual Noncopyable {
96	FamilyFinder *ff;
97	char *pt_server_alignment; // alignment used in pt_server (may differ from 'alignment')
98	int maxRelatives; // max # of relatives to use
99
100	SearchRelativeParams(FamilyFinder ff_, const char pt_server_alignment_, int maxRelatives_)
101	: ff(ff_),
102	pt_server_alignment(strdup(pt_server_alignment_)),
103	maxRelatives(maxRelatives_)
104	{}
105
106	~SearchRelativeParams() {
107	free(pt_server_alignment);
108	delete(ff);
109	}
110
111	FamilyFinder *getFamilyFinder() { return ff; }
112	};
113
114	// --------------------------------------------------------------------------------
115
116	#define GAP_CHAR '-'
117	#define QUALITY_NAME "ASC_ALIGNER_CLIENT_SCORE"
118	#define INSERTS_NAME "AMI_ALIGNER_MASTER_INSERTS"
119
120	#define FA_AWAR_ROOT "faligner/"
121	#define FA_AWAR_TO_ALIGN FA_AWAR_ROOT "what"
122	#define FA_AWAR_REFERENCE FA_AWAR_ROOT "against"
123	#define FA_AWAR_REFERENCE_NAME FA_AWAR_ROOT "sagainst"
124	#define FA_AWAR_RANGE FA_AWAR_ROOT "range"
125	#define FA_AWAR_PROTECTION FA_AWAR_ROOT "protection"
126	#define FA_AWAR_AROUND FA_AWAR_ROOT "around"
127	#define FA_AWAR_MIRROR FA_AWAR_ROOT "mirror"
128	#define FA_AWAR_REPORT FA_AWAR_ROOT "report"
129	#define FA_AWAR_SHOW_GAPS_MESSAGES FA_AWAR_ROOT "show_gaps"
130	#define FA_AWAR_CONTINUE_ON_ERROR FA_AWAR_ROOT "continue_on_error"
131	#define FA_AWAR_ACTION_ON_ERROR FA_AWAR_ROOT "action_on_error"
132	#define FA_AWAR_USE_SECONDARY FA_AWAR_ROOT "use_secondary"
133	#define FA_AWAR_NEXT_RELATIVES FA_AWAR_ROOT "next_relatives"
134	#define FA_AWAR_RELATIVE_RANGE FA_AWAR_ROOT "relrange"
135	#define FA_AWAR_PT_SERVER_ALIGNMENT "tmp/" FA_AWAR_ROOT "relative_ali"
136	#define FA_AWAR_SAI_RANGE_NAME FA_AWAR_ROOT "sai/sainame"
137	#define FA_AWAR_SAI_RANGE_CHARS FA_AWAR_ROOT "sai/chars"
138
139	#define FA_AWAR_ISLAND_HOPPING_ROOT "island_hopping/"
140	#define FA_AWAR_USE_ISLAND_HOPPING FA_AWAR_ISLAND_HOPPING_ROOT "use"
141	#define FA_AWAR_ESTIMATE_BASE_FREQ FA_AWAR_ISLAND_HOPPING_ROOT "estimate_base_freq"
142	#define FA_AWAR_BASE_FREQ_A FA_AWAR_ISLAND_HOPPING_ROOT "base_freq_a"
143	#define FA_AWAR_BASE_FREQ_C FA_AWAR_ISLAND_HOPPING_ROOT "base_freq_c"
144	#define FA_AWAR_BASE_FREQ_G FA_AWAR_ISLAND_HOPPING_ROOT "base_freq_g"
145	#define FA_AWAR_BASE_FREQ_T FA_AWAR_ISLAND_HOPPING_ROOT "base_freq_t"
146	#define FA_AWAR_SUBST_PARA_AC FA_AWAR_ISLAND_HOPPING_ROOT "subst_para_ac"
147	#define FA_AWAR_SUBST_PARA_AG FA_AWAR_ISLAND_HOPPING_ROOT "subst_para_ag"
148	#define FA_AWAR_SUBST_PARA_AT FA_AWAR_ISLAND_HOPPING_ROOT "subst_para_at"
149	#define FA_AWAR_SUBST_PARA_CG FA_AWAR_ISLAND_HOPPING_ROOT "subst_para_cg"
150	#define FA_AWAR_SUBST_PARA_CT FA_AWAR_ISLAND_HOPPING_ROOT "subst_para_ct"
151	#define FA_AWAR_SUBST_PARA_GT FA_AWAR_ISLAND_HOPPING_ROOT "subst_para_gt"
152	#define FA_AWAR_EXPECTED_DISTANCE FA_AWAR_ISLAND_HOPPING_ROOT "expected_dist"
153	#define FA_AWAR_STRUCTURE_SUPPLEMENT FA_AWAR_ISLAND_HOPPING_ROOT "struct_suppl"
154	#define FA_AWAR_THRESHOLD FA_AWAR_ISLAND_HOPPING_ROOT "threshold"
155	#define FA_AWAR_GAP_A FA_AWAR_ISLAND_HOPPING_ROOT "gapa"
156	#define FA_AWAR_GAP_B FA_AWAR_ISLAND_HOPPING_ROOT "gapb"
157	#define FA_AWAR_GAP_C FA_AWAR_ISLAND_HOPPING_ROOT "gapc"
158
159	// --------------------------------------------------------------------------------
160
161	static IslandHopping *island_hopper = NULp; // NULp -> use fast aligner; else use island hopper
162
163	static GB_alignment_type global_alignmentType = GB_AT_UNKNOWN; // type of actually aligned sequence
164
165	static int currentSequenceNumber; // used for counter
166	static int overallSequenceNumber;
167
168	// --------------------------------------------------------------------------------
169
170	inline ARB_ERROR species_not_found(GB_CSTR species_name) {
171	return GBS_global_string("No species '%s' found!", species_name);
172	}
173
174	static ARB_ERROR reverseComplement(GBDATA *gb_species, GB_CSTR ali, int max_protection) {
175	GBDATA *gbd = GBT_find_sequence(gb_species, ali);
176	ARB_ERROR error = NULp;
177
178	if (!gbd) {
179	error = GBS_global_string("No 'data' found for species '%s'", GBT_get_name_or_description(gb_species));
180	}
181	else {
182	int my_protection = GB_read_security_write(gbd);
183
184	if (my_protection<=max_protection) { // ok
185	char *seq = GB_read_string(gbd);
186	int length = GB_read_string_count(gbd);
187	GBDATA *gb_main = GB_get_root(gb_species);
188	GB_alignment_type ali_type = GBT_get_alignment_type(gb_main, ali);
189
190	char T_or_U;
191	error = GBT_determine_T_or_U(ali_type, &T_or_U, "reverse-complement");
192	if (!error) {
193	GBT_reverseComplementNucSequence(seq, length, T_or_U);
194	error = GB_write_string(gbd, seq);
195	}
196	}
197	else { // protection error
198	error = GBS_global_string("Cannot reverse-complement species '%s' because of protection level", GBT_get_name_or_description(gb_species));
199	}
200
201	}
202
203	return error;
204	}
205
206	static void build_reverse_complement(AW_window aw, const AlignDataAccess data_access) {
207	GBDATA *gb_main = data_access->gb_main;
208
209	GB_push_transaction(gb_main);
210
211	AW_root *root = aw->get_root();
212	FA_alignTarget revComplWhat = static_cast<FA_alignTarget>(root->awar(FA_AWAR_TO_ALIGN)->read_int());
213	char *default_alignment = GBT_get_default_alignment(gb_main);
214	GB_CSTR alignment = root->awar_string(AWAR_EDITOR_ALIGNMENT, default_alignment)->read_string();
215	ARB_ERROR error = NULp;
216	int max_protection = root->awar(FA_AWAR_PROTECTION)->read_int();
217
218	switch (revComplWhat) {
219	case FA_CURRENT: { // current species
220	GB_CSTR species_name = root->awar(AWAR_SPECIES_NAME)->read_string();
221	GBDATA *gb_species = GBT_find_species(gb_main, species_name);
222	if (!gb_species) error = species_not_found(species_name);
223	if (!error) error = reverseComplement(gb_species, alignment, max_protection);
224	break;
225	}
226	case FA_MARKED: { // marked species
227	GBDATA *gb_species = GBT_first_marked_species(gb_main);
228
229	if (!gb_species) {
230	error = "There is no marked species";
231	}
232	while (gb_species) {
233	error = reverseComplement(gb_species, alignment, max_protection);
234	if (error) break;
235	gb_species = GBT_next_marked_species(gb_species);
236	}
237	break;
238	}
239	case FA_SELECTED: { // selected species (editor selection!)
240
241	Aligner_get_first_selected_species get_first_selected_species = data_access->get_first_selected_species;
242	Aligner_get_next_selected_species get_next_selected_species = data_access->get_next_selected_species;
243
244	int count = 0;
245	GBDATA *gb_species = get_first_selected_species(&count);
246
247	if (!gb_species) {
248	error = "There is no selected species!";
249	}
250	while (gb_species) {
251	error = reverseComplement(gb_species, alignment, max_protection);
252	if (error) break;
253	gb_species = get_next_selected_species();
254	}
255	break;
256	}
257	default: {
258	fa_assert(0);
259	break;
260	}
261	}
262
263	GB_end_transaction_show_error(gb_main, error, aw_message);
264	}
265
266	// --------------------------------------------------------------------------------
267
268	class AliChange { // describes a local alignment change
269	const CompactedSubSequence& Old;
270	const CompactedSubSequence& New;
271
272	public:
273	AliChange(const CompactedSubSequence& old_, const CompactedSubSequence& new_)
274	: Old(old_), New(new_)
275	{
276	fa_assert(Old.may_refer_to_same_part_as(New));
277	}
278
279	int follow(ExplicitRange& range) const {
280	ExplicitRange compressed(Old.compPosition(range.start()),
281	Old.compPosition(range.end()));
282
283	int exp_start = New.expdPosition(compressed.start());
284	int exp_end = New.expdPosition(compressed.end());
285
286	int gaps_before = New.no_of_gaps_before(compressed.start());
287	int gaps_after = New.no_of_gaps_after(compressed.end());
288
289	fa_assert(gaps_before >= 0);
290	fa_assert(gaps_after >= 0);
291
292	range = ExplicitRange(exp_start-gaps_before, exp_end+gaps_after);
293
294	return compressed.size(); // number of bases
295	}
296	};
297
298	class LooseBases {
299	typedef std::list<ExplicitRange> Ranges;
300
301	Ranges ranges;
302
303	public:
304
305	bool is_empty() const { return ranges.empty(); }
306	void clear() { ranges.clear(); }
307
308	void memorize(ExplicitRange range) {
309	ranges.push_front(range);
310	}
311	ExplicitRange recall() {
312	ExplicitRange range = ranges.front();
313	ranges.pop_front();
314	return range;
315	}
316	int follow_ali_change(const AliChange& change) {
317	// transform positions according to changed alignment (oldSequence -> newSequence)
318	// returns the number of bases contained in 'this'
319	int basecount = 0;
320	if (!is_empty()) {
321	for (Ranges::iterator r = ranges.begin(); r != ranges.end(); ++r) {
322	basecount += change.follow(*r);
323	}
324	}
325	return basecount;
326	}
327	void append(LooseBases& loose) { ranges.splice(ranges.end(), loose.ranges); }
328	int follow_ali_change_and_append(LooseBases& loose, const AliChange& change) {
329	// returns the number of loose bases (added from 'loose')
330	int basecount = loose.follow_ali_change(change);
331	append(loose);
332	return basecount;
333	}
334	};
335
336	static LooseBases unaligned_bases; // if fast_align cannot align (no master bases) -> stores positions here
337
338
339	static const char read_name(GBDATA gbd) {
340	if (!gbd) return "GROUP-CONSENSUS";
341	const char *name = GBT_get_name(gbd);
342	return name ? name : "<unnamed-species>";
343	}
344
345	inline int relatedBases(char base1, char base2) {
346	return baseMatch(base1, base2)==1;
347	}
348
349	inline char alignQuality(char slave, char master) {
350	fa_assert(slave);
351	fa_assert(master);
352
353	char result = '#';
354	if (slave==master) result = '-'; // equal
355	else if (slave==GAP_CHAR) result = '+'; // inserted gap
356	else if (master==GAP_CHAR) result = '+'; // no gap in master
357	else if (relatedBases(slave, master)) result = '~'; // mutation (related bases)
358
359	return result; // mutation (non-related bases)
360	}
361
362	// -------------------------
363	// Debugging stuff
364
365	#ifdef DEBUG
366	static char lstr(const char s, int len) {
367	static int alloc;
368	static char *buffer;
369
370	if (alloc<(len+1)) {
371	if (alloc) free(buffer);
372	ARB_alloc(buffer, alloc=len+100);
373	}
374
375	memcpy(buffer, s, len);
376	buffer[len] = 0;
377
378	return buffer;
379	}
380
381	#define BUFLEN 120
382
383	inline char compareChar(char base1, char base2) {
384	return base1==base2 ? '=' : (relatedBases(base1, base2) ? 'x' : 'X');
385	}
386
387	#if defined(TRACE_COMPRESSED_ALIGNMENT)
388
389	static void dump_n_compare_one(const char seq1, const char seq2, long len, long offset) {
390	fa_assert(len<=BUFLEN);
391	char compare[BUFLEN+1];
392
393	for (long l=0; l<len; l++) {
394	compare[l] = (is_ali_gap(seq1[l]) && is_ali_gap(seq2[l])) ? ' ' : compareChar(seq1[l], seq2[l]);
395	}
396
397	compare[len] = 0;
398
399	printf(" %li '%s'\n", offset, lstr(seq1, len));
400	printf(" %li '%s'\n", offset, lstr(seq2, len));
401	printf(" %li '%s'\n", offset, compare);
402	}
403
404	inline void dump_rest(const char *seq, long len, int idx, long offset) {
405	printf(" Rest von Sequenz %i:\n", idx);
406	while (len>BUFLEN) {
407	printf(" %li '%s'\n", offset, lstr(seq, BUFLEN));
408	seq += BUFLEN;
409	len -= BUFLEN;
410	offset += BUFLEN;
411	}
412
413	fa_assert(len>0);
414	printf(" '%s'\n", lstr(seq, len));
415	}
416
417	static void dump_n_compare(const char text, const char seq1, long len1, const char *seq2, long len2) {
418	long offset = 0;
419
420	printf(" Comparing %s:\n", text);
421
422	while (len1>0 && len2>0) {
423	long done = 0;
424
425	if (len1>=BUFLEN && len2>=BUFLEN) {
426	dump_n_compare_one(seq1, seq2, done=BUFLEN, offset);
427	}
428	else {
429	long min = len1<len2 ? len1 : len2;
430	dump_n_compare_one(seq1, seq2, done=min, offset);
431	}
432
433	seq1 += done;
434	seq2 += done;
435	len1 -= done;
436	len2 -= done;
437	offset += done;
438	}
439
440	if (len1>0) dump_rest(seq1, len1, 1, offset);
441	if (len2>0) dump_rest(seq2, len2, 2, offset);
442	printf(" -------------------\n");
443	}
444
445	static void dump_n_compare(const char *text, const CompactedSubSequence& seq1, const CompactedSubSequence& seq2) {
446	dump_n_compare(text, seq1.text(), seq1.length(), seq2.text(), seq2.length());
447	}
448	#endif // TRACE_COMPRESSED_ALIGNMENT
449
450	#undef BUFLEN
451
452	inline void dumpSeq(const char *seq, long len, long pos) {
453	printf("'%s' ", lstr(seq, len));
454	printf("(Pos=%li,Len=%li)", pos, len);
455	}
456
457	#define dump() \
458	do { \
459	double sig = partSignificance(sequence().length(), slaveSequence.length(), bestLength); \
460	\
461	printf(" Score = %li (Significance=%f)\n" \
462	" Master = ", bestScore, sig); \
463	dumpSeq(bestMasterLeft.text(), bestLength, bestMasterLeft.leftOf()); \
464	printf("\n" \
465	" Slave = "); \
466	dumpSeq(bestSlaveLeft.text(), bestLength, bestSlaveLeft.leftOf()); \
467	printf("\n"); \
468	} while (0)
469
470	#endif //DEBUG
471
472
473	// ------------------------------------------------
474	// INLINE-functions used in fast_align():
475
476	inline double log3(double d) {
477	return log(d)/log(3.0);
478	}
479	inline double partSignificance(long seq1len, long seq2len, long partlen) {
480	// returns log3 of significance of the part
481	// usage: partSignificance(...) < log3(maxAllowedSignificance)
482	return log3((seq1len-partlen)*(seq2len-partlen)) - partlen;
483	}
484
485	inline ARB_ERROR bufferTooSmall() {
486	return "Cannot align - reserved buffer is to small";
487	}
488
489	inline long insertsToNextBase(AlignBuffer *alignBuffer, const SequencePosition& master) {
490	int inserts;
491	int nextBase;
492
493	if (master.rightOf()>0) {
494	nextBase = master.expdPosition();
495	}
496	else {
497	nextBase = master.sequence().expdPosition(master.sequence().length());
498	}
499	inserts = nextBase-alignBuffer->offset();
500
501	return inserts;
502	}
503
504	inline void insertBase(AlignBuffer *alignBuffer,
505	SequencePosition& master, SequencePosition& slave,
506	FastAlignReport *report)
507	{
508	char slaveBase = *slave.text();
509	char masterBase = *master.text();
510
511	alignBuffer->set(slaveBase, alignQuality(slaveBase, masterBase));
512	report->count_aligned_base(slaveBase!=masterBase);
513	++slave;
514	++master;
515	}
516
517	inline void insertSlaveBases(AlignBuffer *alignBuffer,
518	SequencePosition& slave,
519	int length,
520	FastAlignReport *report)
521	{
522	alignBuffer->copy(slave.text(), alignQuality(*slave.text(), GAP_CHAR), length);
523	report->count_unaligned_base(length);
524	slave += length;
525	}
526
527	inline void insertGap(AlignBuffer *alignBuffer,
528	SequencePosition& master,
529	FastAlignReport *report)
530	{
531	char masterBase = *master.text();
532
533	alignBuffer->set(GAP_CHAR, alignQuality(GAP_CHAR, masterBase));
534	report->count_aligned_base(masterBase!=GAP_CHAR);
535	++master;
536	}
537
538	static ARB_ERROR insertClustalValigned(AlignBuffer *alignBuffer,
539	SequencePosition& master,
540	SequencePosition& slave,
541	const char masterAlignment, const char slaveAlignment, long alignmentLength,
542	FastAlignReport *report)
543	{
544	// inserts bases of 'slave' to 'alignBuffer' according to alignment in 'masterAlignment' and 'slaveAlignment'
545	#define ACID '*' // contents of 'masterAlignment' and 'slaveAlignment'
546	#define GAP '-'
547
548	int pos;
549	int baseAtLeft = 0; // TRUE -> last position in alignBuffer contains a base
550
551	for (pos=0; pos<alignmentLength; pos++) {
552	long insert = insertsToNextBase(alignBuffer, master); // in expanded seq
553
554	if (masterAlignment[pos]==ACID) {
555	if (insert>0) {
556	if (insert>alignBuffer->free()) return bufferTooSmall();
557	alignBuffer->set(GAP_CHAR, alignQuality(GAP_CHAR, GAP_CHAR), insert);
558	fa_assert(insertsToNextBase(alignBuffer, master)==0);
559	insert = 0;
560	}
561
562	if (!alignBuffer->free()) return bufferTooSmall();
563	if (slaveAlignment[pos]==ACID) {
564	insertBase(alignBuffer, master, slave, report);
565	baseAtLeft = 1;
566	}
567	else {
568	insertGap(alignBuffer, master, report);
569	baseAtLeft = 0;
570	}
571	}
572	else {
573	int slave_bases;
574
575	fa_assert(masterAlignment[pos]==GAP);
576	for (slave_bases=1; pos+slave_bases<alignmentLength && masterAlignment[pos+slave_bases]==GAP; slave_bases++) {
577	; // count gaps in master (= # of slave bases to insert)
578	}
579	if (!baseAtLeft && insert>slave_bases) {
580	int ins_gaps = insert-slave_bases;
581
582	fa_assert(alignBuffer->free()>=ins_gaps);
583	alignBuffer->set(GAP_CHAR, alignQuality(GAP_CHAR, GAP_CHAR), ins_gaps);
584	insert -= ins_gaps;
585	}
586
587	if (insert<slave_bases) { // master has less gaps than slave bases to insert
588	report->memorize_insertion(master.expdPosition(), slave_bases-insert);
589	}
590	else if (insert>slave_bases) { // master has more gaps than slave bases to insert
591	fa_assert(baseAtLeft);
592	}
593
594	unaligned_bases.memorize(ExplicitRange(slave.expdPosition(), // memorize base positions without counterpart in master
595	slave.expdPosition(slave_bases-1)));
596
597	if (slave_bases>alignBuffer->free()) return bufferTooSmall();
598	insertSlaveBases(alignBuffer, slave, slave_bases, report);
599	pos += slave_bases-1; // -1 compensates for()-increment
600	baseAtLeft = 1;
601	}
602	}
603
604	return NULp;
605
606	#undef GAP
607	#undef ACID
608	}
609
610	static ARB_ERROR insertAligned(AlignBuffer *alignBuffer,
611	SequencePosition& master, SequencePosition& slave, long partLength,
612	FastAlignReport *report)
613	{
614	// insert bases of 'slave' to 'alignBuffer' according to 'master'
615	if (partLength) {
616	long insert = insertsToNextBase(alignBuffer, master);
617
618	fa_assert(partLength>=0);
619
620	if (insert<0) { // insert gaps into master
621	long min_insert = insert;
622
623	report->memorize_insertion(master.expdPosition(), -insert);
624
625	while (insert<0 && partLength) {
626	if (insert<min_insert) min_insert = insert;
627	if (!alignBuffer->free()) {
628	return bufferTooSmall();
629	}
630	insertBase(alignBuffer, master, slave, report);
631	partLength--;
632	insert = insertsToNextBase(alignBuffer, master);
633	}
634
635	fa_assert(partLength>=0);
636	if (partLength==0) { // all inserted
637	return NULp;
638	}
639	}
640
641	fa_assert(insert>=0);
642
643	if (insert>0) { // insert gaps into slave
644	if (insert>alignBuffer->free()) return bufferTooSmall();
645	alignBuffer->set(GAP_CHAR, alignQuality(GAP_CHAR, GAP_CHAR), insert);
646	fa_assert(insertsToNextBase(alignBuffer, master)==0);
647	}
648
649	fa_assert(partLength>=0);
650
651	while (partLength--) {
652	insert = insertsToNextBase(alignBuffer, master);
653
654	fa_assert(insert>=0);
655	if (insert>0) {
656	if (insert>=alignBuffer->free()) return bufferTooSmall();
657	alignBuffer->set(GAP_CHAR, alignQuality(GAP_CHAR, GAP_CHAR), insert);
658	}
659	else {
660	if (!alignBuffer->free()) {
661	return bufferTooSmall();
662	}
663	}
664
665	insertBase(alignBuffer, master, slave, report);
666	}
667	}
668
669	return NULp;
670	}
671
672	static ARB_ERROR cannot_fast_align(const CompactedSubSequence& master, long moffset, long mlength,
673	const CompactedSubSequence& slaveSequence, long soffset, long slength,
674	int max_seq_length,
675	AlignBuffer *alignBuffer,
676	FastAlignReport *report)
677	{
678	const char *mtext = master.text(moffset);
679	const char *stext = slaveSequence.text(soffset);
680	ARB_ERROR error = NULp;
681
682	if (slength) {
683	if (mlength) { // if slave- and master-sequences contain bases, we call the slow aligner
684	#ifdef TRACE_CLUSTAL_DATA
685	printf("ClustalV-Align:\n");
686	printf(" mseq = '%s'\n", lstr(mtext, mlength));
687	printf(" sseq = '%s'\n", lstr(stext, slength));
688	#endif // TRACE_CLUSTAL_DATA
689
690	int is_dna = -1;
691
692	switch (global_alignmentType) {
693	case GB_AT_RNA:
694	case GB_AT_DNA: is_dna = 1; break;
695	case GB_AT_AA: is_dna = 0; break;
696	default: error = "Unknown alignment type - aligner aborted"; break;
697	}
698
699	const char maligned, saligned;
700	int len;
701	if (!error) {
702	int score; // unused
703	error = ClustalV_align(is_dna,
704	1,
705	mtext, mlength, stext, slength,
706	master.gapsBefore(moffset),
707	max_seq_length,
708	maligned, saligned, len, score);
709	}
710
711	if (!error) {
712	#ifdef TRACE_CLUSTAL_DATA
713	printf("ClustalV returns:\n");
714	printf(" maligned = '%s'\n", lstr(maligned, len));
715	printf(" saligned = '%s'\n", lstr(saligned, len));
716	#endif // TRACE_CLUSTAL_DATA
717
718	SequencePosition masterPos(master, moffset);
719	SequencePosition slavePos(slaveSequence, soffset);
720
721	error = insertClustalValigned(alignBuffer, masterPos, slavePos, maligned, saligned, len, report);
722
723	#if (defined(DEBUG) && 0)
724
725	SequencePosition master2(master->sequence(), moffset);
726	SequencePosition slave2(slaveSequence, soffset);
727	char *cmp = new char[len];
728
729	for (int l=0; l<len; l++) {
730	int gaps = 0;
731
732	if (maligned[l]=='*') {
733	maligned[l] = *master2.text();
734	++master2;
735	}
736	else {
737	gaps++;
738	}
739
740	if (saligned[l]=='*') {
741	saligned[l] = *slave2.text();
742	++slave2;
743	}
744	else {
745	gaps++;
746	}
747
748	cmp[l] = gaps \|\| maligned[l]==saligned[l] ? '=' : 'X';
749	}
750
751	printf(" master = '%s'\n", lstr(maligned, len));
752	printf(" slave = '%s'\n", lstr(saligned, len));
753	printf(" '%s'\n", lstr(cmp, len));
754
755	delete [] cmp;
756	#endif
757	}
758	}
759	else { // master is empty here, so we just copy in the slave bases
760	if (slength<=alignBuffer->free()) {
761	unaligned_bases.memorize(ExplicitRange(slaveSequence.expdPosition(soffset),
762	slaveSequence.expdPosition(soffset+slength-1)));
763
764	alignBuffer->copy(slaveSequence.text(soffset), '?', slength); // '?' means not aligned (just inserted)
765	// corrected by at alignBuffer->correctUnalignedPositions()
766	report->count_unaligned_base(slength);
767	}
768	else {
769	error = bufferTooSmall();
770	}
771	}
772	}
773
774	return error;
775	}
776
777	// ------------------------------------
778	// #define's for fast_align()
779
780	#define TEST_BETTER_SCORE() \
781	do { \
782	if (score>bestScore) { \
783	bestScore = score; \
784	bestLength = masterRight.text() - masterLeft.text(); \
785	bestMasterLeft = masterLeft; \
786	bestSlaveLeft = slaveLeft; \
787	} \
788	} while (0)
789
790	#define CAN_SCORE_LEFT() (masterLeft.leftOf() && slaveLeft.leftOf())
791	#define CAN_SCORE_RIGHT() (masterRight.rightOf() && slaveRight.rightOf())
792
793	#define SCORE_LEFT() \
794	do { \
795	score += (--masterLeft).text()==(--slaveLeft).text() ? match : mismatch; \
796	TEST_BETTER_SCORE(); \
797	} while (0)
798
799	#define SCORE_RIGHT() \
800	do { \
801	score += (++masterRight).text()==(++slaveRight).text() ? match : mismatch; \
802	TEST_BETTER_SCORE(); \
803	} while (0)
804
805
806	ARB_ERROR FastSearchSequence::fast_align(const CompactedSubSequence& slaveSequence,
807	AlignBuffer *alignBuffer,
808	int max_seq_length,
809	int match, int mismatch,
810	FastAlignReport *report) const
811	{
812	// aligns 'slaveSequence' to 'this'
813	//
814	// returns
815	// NULp -> all ok -> 'alignBuffer' contains aligned sequence
816	// or error -> failure -> no results
817
818	ARB_ERROR error = NULp;
819	int aligned = 0;
820
821	// set the following #if to zero to test ClustalV-Aligner (not fast_aligner)
822	#if 1
823
824	static double lowSignificance;
825	static int lowSignificanceInitialized;
826
827	if (!lowSignificanceInitialized) {
828	lowSignificance = log3(0.01);
829	lowSignificanceInitialized = 1;
830	}
831
832	SequencePosition slave(slaveSequence);
833	long bestScore=0;
834	SequencePosition bestMasterLeft(sequence());
835	SequencePosition bestSlaveLeft(slaveSequence);
836	long bestLength=0;
837
838	while (slave.rightOf()>=3) { // with all triples of slaveSequence
839	FastSearchOccurrence occurrence(*this, slave.text()); // "search" first
840	SequencePosition rightmostSlave = slave;
841
842	while (occurrence.found()) { // with all occurrences of the triple
843	long score = match*3;
844	SequencePosition masterLeft(occurrence.sequence(), occurrence.offset());
845	SequencePosition masterRight(occurrence.sequence(), occurrence.offset()+3);
846	SequencePosition slaveLeft(slave);
847	SequencePosition slaveRight(slave, 3);
848
849	while (score>0) {
850	if (CAN_SCORE_LEFT()) {
851	SCORE_LEFT();
852	}
853	else {
854	while (score>0 && CAN_SCORE_RIGHT()) {
855	SCORE_RIGHT();
856	}
857	break;
858	}
859
860	if (CAN_SCORE_RIGHT()) {
861	SCORE_RIGHT();
862	}
863	else {
864	while (score>0 && CAN_SCORE_LEFT()) {
865	SCORE_LEFT();
866	}
867	break;
868	}
869	}
870
871	occurrence.gotoNext(); // "search" next
872
873	if (rightmostSlave<slaveRight) {
874	rightmostSlave = slaveRight;
875	rightmostSlave -= 3;
876	}
877	}
878
879	if (rightmostSlave>slave) slave = rightmostSlave;
880	else ++slave;
881	}
882
883	if (bestLength) {
884	double sig = partSignificance(sequence().length(), slaveSequence.length(), bestLength);
885
886	if (sig<lowSignificance) {
887	long masterLeftOf = bestMasterLeft.leftOf();
888	long masterRightStart = masterLeftOf+bestLength;
889	long masterRightOf = bestMasterLeft.rightOf()-bestLength;
890	long slaveLeftOf = bestSlaveLeft.leftOf();
891	long slaveRightStart = slaveLeftOf+bestLength;
892	long slaveRightOf = bestSlaveLeft.rightOf()-bestLength;
893
894	#define MIN_ALIGNMENT_RANGE 4
895
896	if (!error) {
897	if (masterLeftOf >= MIN_ALIGNMENT_RANGE && slaveLeftOf >= MIN_ALIGNMENT_RANGE) {
898	CompactedSubSequence leftCompactedMaster(sequence(), 0, masterLeftOf);
899	FastSearchSequence leftMaster(leftCompactedMaster);
900
901	error = leftMaster.fast_align(CompactedSubSequence(slaveSequence, 0, slaveLeftOf),
902	alignBuffer, max_seq_length, match, mismatch, report);
903	}
904	else if (slaveLeftOf>0) {
905	error = cannot_fast_align(sequence(), 0, masterLeftOf,
906	slaveSequence, 0, slaveLeftOf,
907	max_seq_length, alignBuffer, report);
908	}
909
910	aligned = 1;
911	}
912
913	// align part of slave sequence according to master sequence:
914
915	if (!error) {
916	#if (defined(DEBUG) && 0)
917	long offset = alignBuffer->offset();
918	long used;
919	#endif
920	error = insertAligned(alignBuffer, bestMasterLeft, bestSlaveLeft, bestLength, report);
921	#if (defined(DEBUG) && 0)
922	used = alignBuffer->offset()-offset;
923	printf("aligned '%s' (len=%li, address=%li)\n", lstr(alignBuffer->text()+offset, used), used, long(alignBuffer));
924	#endif
925	aligned = 1;
926	}
927
928	if (!error) {
929	if (masterRightOf >= MIN_ALIGNMENT_RANGE && slaveRightOf >= MIN_ALIGNMENT_RANGE) {
930	CompactedSubSequence rightCompactedMaster(sequence(), masterRightStart, masterRightOf);
931	FastSearchSequence rightMaster(rightCompactedMaster);
932
933	error = rightMaster.fast_align(CompactedSubSequence(slaveSequence, slaveRightStart, slaveRightOf),
934	alignBuffer,
935	max_seq_length, match, mismatch, report);
936	}
937	else if (slaveRightOf>0) {
938	error = cannot_fast_align(sequence(), masterRightStart, masterRightOf,
939	slaveSequence, slaveRightStart, slaveRightOf,
940	max_seq_length, alignBuffer, report);
941	}
942
943	aligned = 1;
944	}
945
946	}
947	}
948
949	#endif
950
951	if (!aligned && !error) {
952	error = cannot_fast_align(sequence(), 0, sequence().length(),
953	slaveSequence, 0, slaveSequence.length(),
954	max_seq_length, alignBuffer, report);
955	}
956
957	return error;
958	}
959
960	#undef dump
961	#undef TEST_BETTER_SCORE
962	#undef CAN_SCORE_LEFT
963	#undef CAN_SCORE_RIGHT
964	#undef SCORE_LEFT
965	#undef SCORE_RIGHT
966
967	inline long calcSequenceChecksum(const char *data, long length) {
968	return GB_checksum(data, length, 1, GAP_CHARS);
969	}
970
971	static CompactedSubSequence readCompactedSequence(GBDATA gb_species,
972	const char *ali,
973	ARB_ERROR *errorPtr,
974	char **dataPtr, // if dataPtr is specified, it will be set to the WHOLE uncompacted sequence
975	long *seqChksum, // may be NULp (of part of sequence)
976	PosRange range) // if !range.is_whole() -> return only part of the sequence
977	{
978	ARB_ERROR error = NULp;
979	GBDATA *gbd;
980	CompactedSubSequence *seq = NULp;
981
982	gbd = GBT_find_sequence(gb_species, ali); // get sequence
983
984	if (gbd) {
985	long length = GB_read_string_count(gbd);
986	char *data = GB_read_string(gbd);
987	long partLength;
988	char *partData;
989
990	if (dataPtr) { // make a copy of the whole uncompacted sequence (returned to caller)
991	*dataPtr = data;
992	}
993
994	int firstColumn = range.start();
995	if (range.is_part()) { // take only part of sequence
996	int lastColumn = range.end();
997
998	fa_assert(firstColumn>=0);
999	fa_assert(firstColumn<=lastColumn);
1000
1001	// include all surrounding gaps (@@@ this might cause problems with partial alignment)
1002	while (firstColumn>0 && is_ali_gap(data[firstColumn-1])) {
1003	firstColumn--;
1004	}
1005	if (lastColumn!=-1) {
1006	while (lastColumn<(length-1) && is_ali_gap(data[lastColumn+1])) lastColumn++;
1007	fa_assert(lastColumn<length);
1008	}
1009
1010	partData = data+firstColumn;
1011	int slen = length-firstColumn;
1012
1013	fa_assert(slen>=0);
1014	fa_assert((size_t)slen==strlen(partData));
1015
1016	if (lastColumn==-1) { // take all till end of sequence
1017	partLength = slen;
1018	}
1019	else {
1020	partLength = lastColumn-firstColumn+1;
1021	if (partLength>slen) partLength = slen; // cut rest, if we have any
1022	}
1023	}
1024	else {
1025	partLength = length;
1026	partData = data;
1027	}
1028
1029	if (!error) {
1030	if (seqChksum) {
1031	*seqChksum = calcSequenceChecksum(partData, partLength);
1032	}
1033
1034	seq = new CompactedSubSequence(partData, partLength, GBT_get_name_or_description(gb_species), firstColumn);
1035	}
1036
1037	if (!dataPtr) { // free sequence only if user has not requested to get it
1038	free(data);
1039	}
1040	}
1041	else {
1042	error = GBS_global_string("No 'data' found for species '%s'", GBT_get_name_or_description(gb_species));
1043	if (dataPtr) *dataPtr = NULp; // (user must not care to free data if we fail)
1044	}
1045
1046	fa_assert(errorPtr);
1047	*errorPtr = error;
1048
1049	return seq;
1050	}
1051
1052	static ARB_ERROR writeStringToAlignment(GBDATA *gb_species, GB_CSTR alignment, GB_CSTR data_name, GB_CSTR str, bool temporary) {
1053	GBDATA *gb_ali = GB_search(gb_species, alignment, GB_DB);
1054	ARB_ERROR error = NULp;
1055	GBDATA *gb_name = GB_search(gb_ali, data_name, GB_STRING);
1056
1057	if (gb_name) {
1058	fa_assert(GB_is_ancestor_of(gb_ali, gb_name));
1059	error = GB_write_string(gb_name, str);
1060	if (temporary && !error) error = GB_set_temporary(gb_name);
1061	}
1062	else {
1063	error = GBS_global_string("Cannot create entry '%s' for '%s'", data_name, GBT_get_name_or_description(gb_species));
1064	}
1065
1066	return error;
1067	}
1068
1069	// --------------------------------------------------------------------------------
1070
1071	static ARB_ERROR alignCompactedTo(CompactedSubSequence *toAlignSequence,
1072	const FastSearchSequence *alignTo,
1073	int max_seq_length,
1074	GB_CSTR alignment,
1075	long toAlignChksum,
1076	GBDATA *gb_toAlign,
1077	GBDATA *gb_alignTo, // may be NULp
1078	const AlignParams& ali_params)
1079	{
1080	// if only part of the sequence should be aligned, then this functions already gets only the part
1081	// (i.o.w.: toAlignSequence, alignTo and toAlignChksum refer to the partial sequence)
1082	AlignBuffer alignBuffer(max_seq_length);
1083	if (ali_params.range.start()>0) {
1084	alignBuffer.set(GAP_CHAR, alignQuality(GAP_CHAR, GAP_CHAR), ali_params.range.start());
1085	}
1086
1087	const char *master_name = read_name(gb_alignTo);
1088
1089	FastAlignReport report(master_name, ali_params.showGapsMessages);
1090
1091	{
1092	static GBDATA *last_gb_toAlign = NULp;
1093	if (gb_toAlign!=last_gb_toAlign) {
1094	last_gb_toAlign = gb_toAlign;
1095	currentSequenceNumber++;
1096	}
1097	}
1098
1099	#ifdef TRACE_COMPRESSED_ALIGNMENT
1100	printf("alignCompactedTo(): master='%s' ", master_name);
1101	printf("slave='%s'\n", toAlignSequence->name());
1102	#endif // TRACE_COMPRESSED_ALIGNMENT
1103
1104	ARB_ERROR error = GB_pop_transaction(gb_toAlign);
1105	if (!error) {
1106	if (island_hopper) {
1107	error = island_hopper->do_align();
1108	if (!error) {
1109	fa_assert(island_hopper->was_aligned());
1110
1111	#ifdef TRACE_ISLANDHOPPER_DATA
1112	printf("Island-Hopper returns:\n");
1113	printf("maligned = '%s'\n", lstr(island_hopper->get_result_ref(), island_hopper->get_result_length()));
1114	printf("saligned = '%s'\n", lstr(island_hopper->get_result(), island_hopper->get_result_length()));
1115	#endif // TRACE_ISLANDHOPPER_DATA
1116
1117	SequencePosition masterPos(alignTo->sequence(), 0);
1118	SequencePosition slavePos(*toAlignSequence, 0);
1119
1120	error = insertClustalValigned(&alignBuffer, masterPos, slavePos, island_hopper->get_result_ref(), island_hopper->get_result(), island_hopper->get_result_length(), &report);
1121	}
1122	}
1123	else {
1124	error = alignTo->fast_align(*toAlignSequence, &alignBuffer, max_seq_length, 2, -10, &report); // <- here we align
1125	}
1126	}
1127
1128	if (!error) {
1129	alignBuffer.correctUnalignedPositions();
1130	if (alignBuffer.free()) {
1131	alignBuffer.set('-', alignQuality(GAP_CHAR, GAP_CHAR), alignBuffer.free()); // rest of alignBuffer is set to '-'
1132	}
1133	alignBuffer.restoreDots(*toAlignSequence);
1134	}
1135
1136	#ifdef TRACE_COMPRESSED_ALIGNMENT
1137	if (!error) {
1138	CompactedSubSequence alignedSlave(alignBuffer.text(), alignBuffer.length(), toAlignSequence->name());
1139	dump_n_compare("reference vs. aligned:", alignTo->sequence(), alignedSlave);
1140	}
1141	#endif // TRACE_COMPRESSED_ALIGNMENT
1142
1143	{
1144	GB_ERROR err = GB_push_transaction(gb_toAlign);
1145	if (!error) error = err;
1146	}
1147
1148	if (!error) {
1149	if (calcSequenceChecksum(alignBuffer.text(), alignBuffer.length())!=toAlignChksum) { // sequence-chksum changed
1150	error = "Internal aligner error (sequence checksum changed) -- aborted";
1151
1152	#ifdef TRACE_COMPRESSED_ALIGNMENT
1153	CompactedSubSequence alignedSlave(alignBuffer.text(), alignBuffer.length(), toAlignSequence->name());
1154	dump_n_compare("Old Slave vs. new Slave", *toAlignSequence, alignedSlave);
1155	#endif // TRACE_COMPRESSED_ALIGNMENT
1156	}
1157	else {
1158	{
1159	GB_topSecurityLevel unsecured(gb_toAlign);
1160	GBDATA *gbd = GBT_add_data(gb_toAlign, alignment, "data", GB_STRING);
1161
1162	if (!gbd) {
1163	error = "Can't find/create sequence data";
1164	}
1165	else {
1166	if (ali_params.range.is_part()) { // we aligned just a part of the sequence
1167	char *buffer = GB_read_string(gbd); // so we read old sequence data
1168	long len = GB_read_string_count(gbd);
1169	if (!buffer) error = GB_await_error();
1170	else {
1171	int lenToCopy = ali_params.range.size();
1172	long wholeChksum = calcSequenceChecksum(buffer, len);
1173
1174	memcpy(buffer+ali_params.range.start(), alignBuffer.text()+ali_params.range.start(), lenToCopy); // copy in the aligned part
1175	// @@@ genau um 1 position zu niedrig
1176
1177	if (calcSequenceChecksum(buffer, len) != wholeChksum) {
1178	error = "Internal aligner error (sequence checksum changed) -- aborted";
1179	# ifdef TRACE_COMPRESSED_ALIGNMENT
1180	char *buffer_org = GB_read_string(gbd);
1181	dump_n_compare("Old seq vs. new seq (slave)", buffer_org, len, buffer, len);
1182	free(buffer_org);
1183	# endif // TRACE_COMPRESSED_ALIGNMENT
1184	}
1185	else {
1186	GB_write_string(gbd, "");
1187	error = GB_write_string(gbd, buffer);
1188	}
1189	}
1190
1191	free(buffer);
1192	}
1193	else {
1194	alignBuffer.setDotsAtEOSequence();
1195	error = GB_write_string(gbd, alignBuffer.text()); // aligned all -> write all
1196	}
1197	}
1198	}
1199
1200	if (!error && ali_params.report != FA_NO_REPORT) {
1201	// create temp-entry for slave containing alignment quality:
1202
1203	error = writeStringToAlignment(gb_toAlign, alignment, QUALITY_NAME, alignBuffer.quality(), ali_params.report==FA_TEMP_REPORT);
1204	if (!error) {
1205	// create temp-entry for master containing insert dots:
1206
1207	int buflen = max_seq_length*2;
1208	char *buffer = ARB_alloc<char>(buflen+1);
1209	char *afterLast = buffer;
1210
1211	if (!buffer) {
1212	error = "out of memory";
1213	}
1214	else {
1215	memset(buffer, '-', buflen);
1216	buffer[buflen] = 0;
1217
1218	const FastAlignInsertion *inserts = report.insertion();
1219	while (inserts) {
1220	memset(buffer+inserts->offset(), '>', inserts->gaps());
1221	afterLast = buffer+inserts->offset()+inserts->gaps();
1222	inserts = inserts->next();
1223	}
1224
1225	*afterLast = 0;
1226	if (gb_alignTo) {
1227	error = writeStringToAlignment(gb_alignTo, alignment, INSERTS_NAME, buffer, ali_params.report==FA_TEMP_REPORT);
1228	}
1229	}
1230	}
1231	}
1232	}
1233	}
1234	return error;
1235	}
1236
1237	ARB_ERROR FastAligner_delete_temp_entries(GBDATA gb_species, const char alignment) {
1238	fa_assert(gb_species);
1239	fa_assert(alignment);
1240
1241	GBDATA *gb_ali = GB_search(gb_species, alignment, GB_FIND);
1242	ARB_ERROR error = NULp;
1243
1244	if (gb_ali) {
1245	GBDATA *gb_name = GB_search(gb_ali, QUALITY_NAME, GB_FIND);
1246	if (gb_name) {
1247	error = GB_delete(gb_name);
1248	#if defined(DEBUG)
1249	printf("deleted '%s' of '%s' (gb_name=%p)\n", QUALITY_NAME, GBT_get_name_or_description(gb_species), gb_name);
1250	#endif
1251	}
1252
1253	if (!error) {
1254	gb_name = GB_search(gb_ali, INSERTS_NAME, GB_FIND);
1255	if (gb_name) {
1256	error = GB_delete(gb_name);
1257	#if defined(DEBUG)
1258	printf("deleted '%s' of '%s' (gb_name=%p)\n", INSERTS_NAME, GBT_get_name_or_description(gb_species), gb_name);
1259	#endif
1260	}
1261	}
1262	}
1263
1264	return error;
1265	}
1266
1267	static ARB_ERROR align_error(ARB_ERROR olderr, GBDATA gb_toAlign, GBDATA gb_alignTo) {
1268	// used by alignTo() and alignToNextRelative() to transform errors coming from subroutines
1269	// can be used by higher functions
1270
1271	const char *name_toAlign = read_name(gb_toAlign);
1272	const char *name_alignTo = read_name(gb_alignTo);
1273
1274	fa_assert(olderr);
1275
1276	return GBS_global_string("Error while aligning '%s' to '%s':\n%s",
1277	name_toAlign, name_alignTo, olderr.deliver());
1278	}
1279
1280	static ARB_ERROR alignTo(GBDATA *gb_toAlign,
1281	GB_CSTR alignment,
1282	const FastSearchSequence *alignTo,
1283	GBDATA *gb_alignTo, // may be NULp
1284	int max_seq_length,
1285	const AlignParams& ali_params)
1286	{
1287	ARB_ERROR error = NULp;
1288	long chksum = -1;
1289
1290	CompactedSubSequence *toAlignSequence = readCompactedSequence(gb_toAlign, alignment, &error, NULp, &chksum, ali_params.range);
1291
1292	if (island_hopper) {
1293	GBDATA *gb_seq = GBT_find_sequence(gb_toAlign, alignment); // get sequence
1294	if (gb_seq) {
1295	long length = GB_read_string_count(gb_seq);
1296	const char *data = GB_read_char_pntr(gb_seq);
1297
1298	island_hopper->set_toAlign_sequence(data);
1299	island_hopper->set_alignment_length(length);
1300	}
1301	}
1302
1303
1304
1305	if (!error) {
1306	error = alignCompactedTo(toAlignSequence, alignTo, max_seq_length, alignment, chksum, gb_toAlign, gb_alignTo, ali_params);
1307	if (error) error = align_error(error, gb_toAlign, gb_alignTo);
1308	delete toAlignSequence;
1309	}
1310
1311	return error;
1312	}
1313
1314	static ARB_ERROR alignToGroupConsensus(GBDATA *gb_toAlign,
1315	GB_CSTR alignment,
1316	Aligner_get_consensus_func get_consensus,
1317	int max_seq_length,
1318	const AlignParams& ali_params)
1319	{
1320	ARB_ERROR error = NULp;
1321	char *consensus = get_consensus(read_name(gb_toAlign), ali_params.range);
1322	size_t cons_len = strlen(consensus);
1323
1324	fa_assert(cons_len);
1325
1326	for (size_t i = 0; i<cons_len; ++i) { // translate consensus to be accepted by aligner
1327	switch (consensus[i]) {
1328	case '=': consensus[i] = '-'; break;
1329	default: break;
1330	}
1331	}
1332
1333	CompactedSubSequence compacted(consensus, cons_len, "group consensus", ali_params.range.start());
1334
1335	{
1336	FastSearchSequence fast(compacted);
1337	error = alignTo(gb_toAlign, alignment, &fast, NULp, max_seq_length, ali_params);
1338	}
1339
1340	free(consensus);
1341
1342	return error;
1343	}
1344
1345	static void appendNameAndUsedBasePositions(char *toString, GBDATA gb_species, int usedBasePositions) {
1346	// if usedBasePositions == -1 -> unknown;
1347
1348	char *currInfo;
1349	if (usedBasePositions<0) {
1350	currInfo = strdup(GBT_get_name_or_description(gb_species));
1351	}
1352	else {
1353	fa_assert(usedBasePositions>0); // otherwise it should NOT be announced here!
1354	currInfo = GBS_global_string_copy("%s:%i", GBT_get_name_or_description(gb_species), usedBasePositions);
1355	}
1356
1357	char *newString = NULp;
1358	if (*toString) {
1359	newString = GBS_global_string_copy("%s, %s", *toString, currInfo);
1360	}
1361	else {
1362	newString = currInfo;
1363	currInfo = NULp; // don't free
1364	}
1365
1366	freeset(*toString, newString);
1367	free(currInfo);
1368	}
1369
1370	inline int min(int i, int j) { return i<j ? i : j; }
1371
1372	static ARB_ERROR alignToNextRelative(SearchRelativeParams& relSearch,
1373	int max_seq_length,
1374	FA_turn turnAllowed,
1375	GB_CSTR alignment,
1376	GBDATA *gb_toAlign,
1377	const AlignParams& ali_params)
1378	{
1379	CompactedSubSequence *toAlignSequence = NULp;
1380	bool use_different_pt_server_alignment = 0 != strcmp(relSearch.pt_server_alignment, alignment);
1381
1382	ARB_ERROR error;
1383	long chksum;
1384	int relativesToTest = relSearch.maxRelatives*2; // get more relatives from pt-server (needed when use_different_pt_server_alignment == true)
1385	char *nearestRelative = new char[relativesToTest+1];
1386	int next_relatives;
1387	int i;
1388	GBDATA *gb_main = GB_get_root(gb_toAlign);
1389
1390	if (use_different_pt_server_alignment) {
1391	turnAllowed = FA_TURN_NEVER; // makes no sense if we're using a different alignment for the pt_server
1392	}
1393
1394	for (next_relatives=0; next_relatives<relativesToTest; next_relatives++) {
1395	nearestRelative[next_relatives] = NULp;
1396	}
1397	next_relatives = 0;
1398
1399	bool restart = true;
1400	while (restart) {
1401	restart = false;
1402
1403	char *findRelsBySeq = NULp;
1404	if (use_different_pt_server_alignment) {
1405	toAlignSequence = readCompactedSequence(gb_toAlign, alignment, &error, NULp, &chksum, ali_params.range);
1406
1407	GBDATA *gbd = GBT_find_sequence(gb_toAlign, relSearch.pt_server_alignment); // use a different alignment for next relative search
1408	if (!gbd) {
1409	error = GBS_global_string("Species '%s' has no data in alignment '%s'", GBT_get_name_or_description(gb_toAlign), relSearch.pt_server_alignment);
1410	}
1411	else {
1412	findRelsBySeq = GB_read_string(gbd);
1413	}
1414	}
1415	else {
1416	toAlignSequence = readCompactedSequence(gb_toAlign, alignment, &error, &findRelsBySeq, &chksum, ali_params.range);
1417	}
1418
1419	if (error) {
1420	delete toAlignSequence;
1421	return error; // @@@ leaks ?
1422	}
1423
1424	while (next_relatives) {
1425	next_relatives--;
1426	freenull(nearestRelative[next_relatives]);
1427	}
1428
1429	{
1430	// find relatives
1431	FamilyFinder *familyFinder = relSearch.getFamilyFinder();
1432	const PosRange& range = familyFinder->get_TargetRange();
1433
1434	if (range.is_part()) {
1435	range.copy_corresponding_part(findRelsBySeq, findRelsBySeq, strlen(findRelsBySeq));
1436	turnAllowed = FA_TURN_NEVER; // makes no sense if we're using partial relative search
1437	}
1438
1439	error = familyFinder->searchFamily(findRelsBySeq, FF_FORWARD, relativesToTest+1, 0); // @@@ make min_score configurable
1440
1441	// @@@ case where no relative found (due to min score) handle how ? abort ? warn ?
1442
1443	double bestScore = 0;
1444	if (!error) {
1445	#if defined(DEBUG)
1446	double lastScore = -1;
1447	#if defined(TRACE_RELATIVES)
1448	fprintf(stderr, "List of relatives used for '%s':\n", GBT_get_name_or_description(gb_toAlign));
1449	#endif // TRACE_RELATIVES
1450	#endif // DEBUG
1451	for (const FamilyList *fl = familyFinder->getFamilyList(); fl; fl=fl->next) {
1452	if (strcmp(toAlignSequence->name(), fl->name)!=0) {
1453	if (GBT_find_species(gb_main, fl->name)) { // @@@
1454	double thisScore = familyFinder->uses_rel_matches() ? fl->rel_matches : fl->matches;
1455	#if defined(DEBUG)
1456	// check whether family list is sorted correctly
1457	fa_assert(lastScore < 0 \|\| lastScore >= thisScore);
1458	lastScore = thisScore;
1459	#if defined(TRACE_RELATIVES)
1460	fprintf(stderr, "- %s (%5.2f)\n", fl->name, thisScore);
1461	#endif // TRACE_RELATIVES
1462	#endif // DEBUG
1463	if (thisScore>=bestScore) bestScore = thisScore;
1464	if (next_relatives<(relativesToTest+1)) {
1465	nearestRelative[next_relatives] = strdup(fl->name);
1466	next_relatives++;
1467	}
1468	}
1469	}
1470	}
1471	}
1472
1473	if (!error && turnAllowed != FA_TURN_NEVER) { // test if mirrored sequence has better relatives
1474	char *mirroredSequence = strdup(findRelsBySeq);
1475	long length = strlen(mirroredSequence);
1476	double bestMirroredScore = 0;
1477
1478	char T_or_U;
1479	error = GBT_determine_T_or_U(global_alignmentType, &T_or_U, "reverse-complement");
1480	if (!error) {
1481	GBT_reverseComplementNucSequence(mirroredSequence, length, T_or_U);
1482	error = familyFinder->searchFamily(mirroredSequence, FF_FORWARD, relativesToTest+1, 0); // @@@ make min_score configurable
1483	}
1484	if (!error) {
1485	#if defined(DEBUG)
1486	double lastScore = -1;
1487	#if defined(TRACE_RELATIVES)
1488	fprintf(stderr, "List of relatives used for '%s' (turned seq):\n", GBT_get_name_or_description(gb_toAlign));
1489	#endif // TRACE_RELATIVES
1490	#endif // DEBUG
1491	for (const FamilyList *fl = familyFinder->getFamilyList(); fl; fl = fl->next) {
1492	double thisScore = familyFinder->uses_rel_matches() ? fl->rel_matches : fl->matches;
1493	#if defined(DEBUG)
1494	// check whether family list is sorted correctly
1495	fa_assert(lastScore < 0 \|\| lastScore >= thisScore);
1496	lastScore = thisScore;
1497	#if defined(TRACE_RELATIVES)
1498	fprintf(stderr, "- %s (%5.2f)\n", fl->name, thisScore);
1499	#endif // TRACE_RELATIVES
1500	#endif // DEBUG
1501	if (thisScore >= bestMirroredScore) {
1502	if (strcmp(toAlignSequence->name(), fl->name)!=0) {
1503	if (GBT_find_species(gb_main, fl->name)) bestMirroredScore = thisScore; // @@@
1504	}
1505	}
1506	}
1507	}
1508
1509	int turnIt = 0;
1510	if (bestMirroredScore>bestScore) {
1511	if (turnAllowed==FA_TURN_INTERACTIVE) {
1512	const char *message;
1513	if (familyFinder->uses_rel_matches()) {
1514	message = GBS_global_string("'%s' seems to be the other way round (score: %.1f%%, score if turned: %.1f%%)",
1515	toAlignSequence->name(), bestScore100, bestMirroredScore100);
1516	}
1517	else {
1518	message = GBS_global_string("'%s' seems to be the other way round (score: %li, score if turned: %li)",
1519	toAlignSequence->name(), long(bestScore+.5), long(bestMirroredScore+.5));
1520	}
1521	turnIt = aw_question("fastali_turn_sequence", message, "Turn sequence,Leave sequence alone")==0;
1522	}
1523	else {
1524	fa_assert(turnAllowed == FA_TURN_ALWAYS);
1525	turnIt = 1;
1526	#if defined(TRACE_RELATIVES)
1527	fprintf(stderr, "Using turned sequence!\n");
1528	#endif // TRACE_RELATIVES
1529	}
1530	}
1531
1532	if (turnIt) { // write mirrored sequence
1533	GBDATA *gbd = GBT_find_sequence(gb_toAlign, alignment);
1534	GB_topSecurityLevel unsecured(gbd);
1535	error = GB_write_string(gbd, mirroredSequence);
1536	if (!error) {
1537	delete toAlignSequence;
1538	restart = true; // continue while loop
1539	}
1540	}
1541
1542	free(mirroredSequence);
1543	}
1544	}
1545	free(findRelsBySeq);
1546	}
1547
1548	if (!error) {
1549	if (!next_relatives) {
1550	char warning[200];
1551	sprintf(warning, "No relative found for '%s'", toAlignSequence->name());
1552	aw_message(warning);
1553	}
1554	else {
1555	// assuming relatives are sorted! (nearest to farthest)
1556
1557	// get data pointers
1558	typedef GBDATA *GBDATAP;
1559	GBDATAP *gb_reference = new GBDATAP[relSearch.maxRelatives];
1560
1561	{
1562	for (i=0; i<relSearch.maxRelatives && i<next_relatives; i++) {
1563	GBDATA *gb_species = GBT_find_species(gb_main, nearestRelative[i]);
1564	if (!gb_species) { // pt-server seems not up to date!
1565	error = species_not_found(nearestRelative[i]);
1566	break;
1567	}
1568
1569	GBDATA *gb_sequence = GBT_find_sequence(gb_species, alignment);
1570	if (gb_sequence) { // if relative has sequence data in the current alignment ..
1571	gb_reference[i] = gb_species;
1572	}
1573	else { // remove this relative
1574	free(nearestRelative[i]);
1575	for (int j = i+1; j<next_relatives; ++j) {
1576	nearestRelative[j-1] = nearestRelative[j];
1577	}
1578	next_relatives--;
1579	nearestRelative[next_relatives] = NULp;
1580	i--; // redo same index
1581	}
1582	}
1583
1584	// delete superfluous relatives
1585	for (; i<next_relatives; ++i) freenull(nearestRelative[i]);
1586
1587	if (next_relatives>relSearch.maxRelatives) {
1588	next_relatives = relSearch.maxRelatives;
1589	}
1590	}
1591
1592	// align
1593
1594	if (!error) {
1595	CompactedSubSequence *alignToSequence = readCompactedSequence(gb_reference[0], alignment, &error, NULp, NULp, ali_params.range);
1596	fa_assert(alignToSequence);
1597
1598	if (island_hopper) {
1599	GBDATA *gb_ref = GBT_find_sequence(gb_reference[0], alignment); // get reference sequence
1600	GBDATA *gb_align = GBT_find_sequence(gb_toAlign, alignment); // get sequence to align
1601
1602	if (gb_ref && gb_align) {
1603	long length_ref = GB_read_string_count(gb_ref);
1604	const char *data;
1605
1606	data = GB_read_char_pntr(gb_ref);
1607	island_hopper->set_ref_sequence(data);
1608
1609	data = GB_read_char_pntr(gb_align);
1610	island_hopper->set_toAlign_sequence(data);
1611
1612	island_hopper->set_alignment_length(length_ref);
1613	}
1614	}
1615
1616	{
1617	FastSearchSequence referenceFastSeq(*alignToSequence);
1618
1619	error = alignCompactedTo(toAlignSequence, &referenceFastSeq,
1620	max_seq_length, alignment, chksum,
1621	gb_toAlign, gb_reference[0], ali_params);
1622	}
1623
1624	if (error) {
1625	error = align_error(error, gb_toAlign, gb_reference[0]);
1626	}
1627	else {
1628	char *used_relatives = NULp;
1629
1630	if (unaligned_bases.is_empty()) {
1631	appendNameAndUsedBasePositions(&used_relatives, gb_reference[0], -1);
1632	}
1633	else {
1634	if (island_hopper) {
1635	appendNameAndUsedBasePositions(&used_relatives, gb_reference[0], -1);
1636	if (next_relatives>1) error = "Island hopping uses only one relative";
1637	}
1638	else {
1639	LooseBases loose;
1640	LooseBases loose_for_next_relative;
1641
1642	int unaligned_positions;
1643	{
1644	CompactedSubSequence *alignedSequence = readCompactedSequence(gb_toAlign, alignment, &error, NULp, NULp, ali_params.range);
1645
1646	unaligned_positions = loose.follow_ali_change_and_append(unaligned_bases, AliChange(toAlignSequence, alignedSequence));
1647	// now loose holds the unaligned (and recalculated) parts from last relative
1648	delete alignedSequence;
1649	}
1650
1651	if (!error) {
1652	int toalign_positions = toAlignSequence->length();
1653	if (unaligned_positions<toalign_positions) {
1654	appendNameAndUsedBasePositions(&used_relatives, gb_reference[0], toalign_positions-unaligned_positions);
1655	}
1656	}
1657
1658	for (i=1; i<next_relatives && !error; i++) {
1659	loose.append(loose_for_next_relative);
1660	int unaligned_positions_for_next = 0;
1661
1662	while (!loose.is_empty() && !error) {
1663	ExplicitRange partRange(intersection(loose.recall(), ali_params.range));
1664	CompactedSubSequence *alignToPart = readCompactedSequence(gb_reference[i], alignment, &error, NULp, NULp, partRange);
1665
1666	if (!error) {
1667	long part_chksum;
1668	CompactedSubSequence *toAlignPart = readCompactedSequence(gb_toAlign, alignment, &error, NULp, &part_chksum, partRange);
1669
1670	fa_assert(contradicted(error, toAlignPart));
1671
1672	if (!error) {
1673	AlignParams loose_ali_params = { ali_params.report, ali_params.showGapsMessages, partRange };
1674
1675	FastSearchSequence referenceFastSeq(*alignToPart);
1676	error = alignCompactedTo(toAlignPart, &referenceFastSeq,
1677	max_seq_length, alignment, part_chksum,
1678	gb_toAlign, gb_reference[i], loose_ali_params);
1679	if (!error) {
1680	CompactedSubSequence *alignedPart = readCompactedSequence(gb_toAlign, alignment, &error, NULp, NULp, partRange);
1681	if (!error) {
1682	unaligned_positions_for_next += loose_for_next_relative.follow_ali_change_and_append(unaligned_bases, AliChange(toAlignPart, alignedPart));
1683	}
1684	delete alignedPart;
1685	}
1686	}
1687	delete toAlignPart;
1688	}
1689	delete alignToPart;
1690	}
1691
1692	if (!error) {
1693	fa_assert(unaligned_positions_for_next <= unaligned_positions); // means: number of unaligned positions has increased by use of relative
1694	if (unaligned_positions_for_next<unaligned_positions) {
1695	appendNameAndUsedBasePositions(&used_relatives, gb_reference[i], unaligned_positions-unaligned_positions_for_next);
1696	unaligned_positions = unaligned_positions_for_next;
1697	}
1698	}
1699	}
1700	}
1701	}
1702
1703	if (!error) { // write used relatives to db-field
1704	error = GBT_write_string(gb_toAlign, "used_rels", used_relatives);
1705	}
1706	free(used_relatives);
1707	}
1708
1709	delete alignToSequence;
1710	}
1711
1712	delete [] gb_reference;
1713	}
1714	}
1715
1716	delete toAlignSequence;
1717
1718	for (i=0; i<next_relatives; i++) freenull(nearestRelative[i]);
1719	delete [] nearestRelative;
1720
1721	return error;
1722	}
1723
1724	// ------------------------
1725	// AlignmentReference
1726
1727	class AlignmentReference : virtual Noncopyable {
1728	GB_CSTR alignment;
1729	int max_seq_length;
1730	const AlignParams& ali_params;
1731
1732	public:
1733	AlignmentReference(GB_CSTR alignment_,
1734	int max_seq_length_,
1735	const AlignParams& ali_params_)
1736	: alignment(alignment_),
1737	max_seq_length(max_seq_length_),
1738	ali_params(ali_params_)
1739	{}
1740	virtual ~AlignmentReference() {}
1741
1742	virtual ARB_ERROR align_to(GBDATA *gb_toalign) const = 0;
1743
1744	GB_CSTR get_alignment() const { return alignment; }
1745	int get_max_seq_length() const { return max_seq_length; }
1746	const AlignParams& get_ali_params() const { return ali_params; }
1747	};
1748
1749
1750	// @@@ make alignTo a member of ExplicitReference (or of AlignmentReference)
1751	// @@@ let alignToGroupConsensus and alignToNextRelative use ExplicitReference
1752
1753
1754	class ExplicitReference: public AlignmentReference { // derived from a Noncopyable
1755	const FastSearchSequence *targetSequence;
1756	GBDATA *gb_alignTo;
1757
1758	public:
1759	ExplicitReference(GB_CSTR alignment_,
1760	const FastSearchSequence *targetSequence_,
1761	GBDATA *gb_alignTo_,
1762	int max_seq_length_,
1763	const AlignParams& ali_params_)
1764	: AlignmentReference(alignment_, max_seq_length_, ali_params_),
1765	targetSequence(targetSequence_),
1766	gb_alignTo(gb_alignTo_)
1767	{}
1768
1769	ARB_ERROR align_to(GBDATA *gb_toalign) const OVERRIDE {
1770	return alignTo(gb_toalign, get_alignment(), targetSequence, gb_alignTo, get_max_seq_length(), get_ali_params());
1771	}
1772	};
1773
1774	// @@@ make alignToGroupConsensus a member of ConsensusReference
1775
1776	class ConsensusReference: public AlignmentReference {
1777	Aligner_get_consensus_func get_consensus;
1778
1779	public:
1780	ConsensusReference(GB_CSTR alignment_,
1781	Aligner_get_consensus_func get_consensus_,
1782	int max_seq_length_,
1783	const AlignParams& ali_params_)
1784	: AlignmentReference(alignment_, max_seq_length_, ali_params_),
1785	get_consensus(get_consensus_)
1786	{}
1787
1788	ARB_ERROR align_to(GBDATA *gb_toalign) const OVERRIDE {
1789	return alignToGroupConsensus(gb_toalign, get_alignment(), get_consensus, get_max_seq_length(), get_ali_params());
1790	}
1791	};
1792
1793	// @@@ make alignToNextRelative a member of SearchRelativesReference
1794
1795	class SearchRelativesReference: public AlignmentReference {
1796	SearchRelativeParams& relSearch;
1797	FA_turn turnAllowed;
1798
1799	public:
1800	SearchRelativesReference(SearchRelativeParams& relSearch_,
1801	int max_seq_length_,
1802	FA_turn turnAllowed_,
1803	GB_CSTR alignment_,
1804	const AlignParams& ali_params_)
1805	: AlignmentReference(alignment_, max_seq_length_, ali_params_),
1806	relSearch(relSearch_),
1807	turnAllowed(turnAllowed_)
1808	{}
1809
1810	ARB_ERROR align_to(GBDATA *gb_toalign) const OVERRIDE {
1811	return alignToNextRelative(relSearch, get_max_seq_length(), turnAllowed, get_alignment(), gb_toalign, get_ali_params());
1812	}
1813	};
1814
1815
1816	// ----------------
1817	// Aligner
1818
1819	class Aligner : virtual Noncopyable {
1820	GBDATA *gb_main;
1821
1822	// define alignment target(s):
1823	FA_alignTarget alignWhat;
1824	GB_CSTR alignment; // name of alignment to use (if NULp -> use default)
1825	GB_CSTR toalign; // name of species to align (used if alignWhat == FA_CURRENT)
1826	Aligner_get_first_selected_species get_first_selected_species; // used if alignWhat == FA_SELECTED
1827	Aligner_get_next_selected_species get_next_selected_species; // --- "" ---
1828
1829	// define reference sequence(s):
1830	GB_CSTR reference; // name of reference species
1831	Aligner_get_consensus_func get_consensus; // function to get consensus seq
1832	SearchRelativeParams& relSearch; // params to search for relatives
1833
1834	// general params:
1835	FA_turn turnAllowed;
1836	const AlignParams& ali_params;
1837	int maxProtection; // protection level
1838
1839	// -------------------- new members
1840	int wasNotAllowedToAlign; // number of failures caused by wrong protection
1841	int err_count; // count errors
1842	bool continue_on_error; /* true -> run single alignments in separate transactions.
1843	* If one target fails, continue with rest.
1844	* false -> run all in one transaction
1845	* One fails -> all fail!
1846	*/
1847	FA_errorAction error_action;
1848
1849	typedef std::list<GBDATA*> GBDATAlist;
1850	GBDATAlist species_to_mark; // species that will be marked after aligning
1851
1852	ARB_ERROR alignToReference(GBDATA *gb_toalign, const AlignmentReference& ref);
1853	ARB_ERROR alignTargetsToReference(const AlignmentReference& ref, GBDATA *gb_species_data);
1854
1855	ARB_ERROR alignToExplicitReference(GBDATA *gb_species_data, int max_seq_length);
1856	ARB_ERROR alignToConsensus(GBDATA *gb_species_data, int max_seq_length);
1857	ARB_ERROR alignToRelatives(GBDATA *gb_species_data, int max_seq_length);
1858
1859	void triggerAction(GBDATA *gb_species, bool has_been_aligned) {
1860	bool mark = false;
1861	switch (error_action) {
1862	case FA_MARK_FAILED: mark = !has_been_aligned; break;
1863	case FA_MARK_ALIGNED: mark = has_been_aligned; break;
1864	case FA_NO_ACTION: mark = false; break;
1865	}
1866	if (mark) species_to_mark.push_back(gb_species);
1867	}
1868
1869	public:
1870
1871	// @@@ pass AlignmentReference from caller (replacing reference parameters)
1872
1873	Aligner(GBDATA *gb_main_,
1874
1875	// define alignment target(s):
1876	FA_alignTarget alignWhat_,
1877	GB_CSTR alignment_, // name of alignment to use (if NULp -> use default)
1878	GB_CSTR toalign_, // name of species to align (used if alignWhat == FA_CURRENT)
1879	Aligner_get_first_selected_species get_first_selected_species_, // used if alignWhat == FA_SELECTED
1880	Aligner_get_next_selected_species get_next_selected_species_, // --- "" ---
1881
1882	// define reference sequence(s):
1883	GB_CSTR reference_, // name of reference species
1884	Aligner_get_consensus_func get_consensus_, // function to get consensus seq
1885	SearchRelativeParams& relSearch_, // params to search for relatives
1886
1887	// general params:
1888	FA_turn turnAllowed_,
1889	const AlignParams& ali_params_,
1890	int maxProtection_, // protection level
1891	bool continue_on_error_,
1892	FA_errorAction error_action_)
1893	: gb_main(gb_main_),
1894	alignWhat(alignWhat_),
1895	alignment(alignment_),
1896	toalign(toalign_),
1897	get_first_selected_species(get_first_selected_species_),
1898	get_next_selected_species(get_next_selected_species_),
1899	reference(reference_),
1900	get_consensus(get_consensus_),
1901	relSearch(relSearch_),
1902	turnAllowed(turnAllowed_),
1903	ali_params(ali_params_),
1904	maxProtection(maxProtection_),
1905	wasNotAllowedToAlign(0),
1906	err_count(0),
1907	continue_on_error(continue_on_error_),
1908	error_action(continue_on_error ? error_action_ : FA_NO_ACTION)
1909	{}
1910
1911	ARB_ERROR run();
1912	};
1913
1914	ARB_ERROR Aligner::alignToReference(GBDATA *gb_toalign, const AlignmentReference& ref) {
1915	int myProtection = GB_read_security_write(GBT_find_sequence(gb_toalign, alignment));
1916	ARB_ERROR error;
1917
1918	if (myProtection<=maxProtection) {
1919	// Depending on 'continue_on_error' we either
1920	// * stop aligning if an error occurs or
1921	// * run the alignment of each species in its own transaction, ignore but report errors and
1922	// optionally mark aligned or failed species.
1923
1924	if (continue_on_error) {
1925	fa_assert(GB_get_transaction_level(gb_main) == 1);
1926	error = GB_end_transaction(gb_main, error); // end global transaction
1927	}
1928
1929	if (!error) {
1930	error = GB_push_transaction(gb_main);
1931	if (!error) error = ref.align_to(gb_toalign);
1932	error = GB_end_transaction(gb_main, error);
1933
1934	if (error) err_count++;
1935	triggerAction(gb_toalign, !error);
1936	}
1937
1938	if (continue_on_error) {
1939	if (error) {
1940	GB_warning(error.deliver());
1941	error = NULp;
1942	}
1943	error = GB_begin_transaction(gb_main); // re-open global transaction
1944	}
1945	}
1946	else {
1947	wasNotAllowedToAlign++;
1948	triggerAction(gb_toalign, false);
1949	}
1950
1951	return error;
1952	}
1953
1954	ARB_ERROR Aligner::alignTargetsToReference(const AlignmentReference& ref, GBDATA *gb_species_data) {
1955	ARB_ERROR error;
1956
1957	fa_assert(GB_get_transaction_level(gb_main) == 1);
1958
1959	switch (alignWhat) {
1960	case FA_CURRENT: { // align one sequence
1961	fa_assert(toalign);
1962
1963	GBDATA *gb_toalign = GBT_find_species_rel_species_data(gb_species_data, toalign);
1964	if (!gb_toalign) {
1965	error = species_not_found(toalign);
1966	}
1967	else {
1968	currentSequenceNumber = overallSequenceNumber = 1;
1969	error = alignToReference(gb_toalign, ref);
1970	}
1971	break;
1972	}
1973	case FA_MARKED: { // align all marked sequences
1974	int count = GBT_count_marked_species(gb_main);
1975	GBDATA *gb_species = GBT_first_marked_species_rel_species_data(gb_species_data);
1976
1977	arb_progress progress("Aligning marked species", long(count));
1978	progress.auto_subtitles("Species");
1979
1980	currentSequenceNumber = 1;
1981	overallSequenceNumber = count;
1982
1983	while (gb_species && !error) {
1984	error = alignToReference(gb_species, ref);
1985	progress.inc_and_check_user_abort(error);
1986	gb_species = GBT_next_marked_species(gb_species);
1987	}
1988	break;
1989	}
1990	case FA_SELECTED: { // align all selected species
1991	int count;
1992	GBDATA *gb_species = get_first_selected_species(&count);
1993
1994
1995	currentSequenceNumber = 1;
1996	overallSequenceNumber = count;
1997
1998	if (!gb_species) {
1999	aw_message("There is no selected species!");
2000	}
2001	else {
2002	arb_progress progress("Aligning selected species", long(count));
2003	progress.auto_subtitles("Species");
2004
2005	while (gb_species && !error) {
2006	error = alignToReference(gb_species, ref);
2007	progress.inc_and_check_user_abort(error);
2008	gb_species = get_next_selected_species();
2009	}
2010	}
2011	break;
2012	}
2013	}
2014
2015	fa_assert(GB_get_transaction_level(gb_main) == 1);
2016	return error;
2017	}
2018
2019	ARB_ERROR Aligner::alignToExplicitReference(GBDATA *gb_species_data, int max_seq_length) {
2020	ARB_ERROR error;
2021	GBDATA *gb_reference = GBT_find_species_rel_species_data(gb_species_data, reference);
2022
2023	if (!gb_reference) {
2024	error = species_not_found(reference);
2025	}
2026	else {
2027	long referenceChksum;
2028	CompactedSubSequence *referenceSeq = readCompactedSequence(gb_reference, alignment, &error, NULp, &referenceChksum, ali_params.range);
2029
2030	if (island_hopper) {
2031	// @@@ setting island_hopper reference has to be done in called function (seems that it is NOT done for alignToConsensus and alignToRelatives). First get tests in place!
2032
2033	GBDATA *gb_seq = GBT_find_sequence(gb_reference, alignment); // get sequence
2034	if (gb_seq) {
2035	long length = GB_read_string_count(gb_seq);
2036	const char *data = GB_read_char_pntr(gb_seq);
2037
2038	island_hopper->set_ref_sequence(data);
2039	island_hopper->set_alignment_length(length);
2040	}
2041	}
2042
2043
2044	if (!error) {
2045	// @@@ do not pass FastSearchSequence to ExplicitReference, instead pass sequence and length (ExplicitReference shall create it itself)
2046
2047	FastSearchSequence referenceFastSeq(*referenceSeq);
2048	ExplicitReference target(alignment, &referenceFastSeq, gb_reference, max_seq_length, ali_params);
2049
2050	error = alignTargetsToReference(target, gb_species_data);
2051	}
2052	delete referenceSeq;
2053	}
2054	return error;
2055	}
2056
2057	ARB_ERROR Aligner::alignToConsensus(GBDATA *gb_species_data, int max_seq_length) {
2058	return alignTargetsToReference(ConsensusReference(alignment, get_consensus, max_seq_length, ali_params),
2059	gb_species_data);
2060	}
2061
2062	ARB_ERROR Aligner::alignToRelatives(GBDATA *gb_species_data, int max_seq_length) {
2063	return alignTargetsToReference(SearchRelativesReference(relSearch, max_seq_length, turnAllowed, alignment, ali_params),
2064	gb_species_data);
2065	}
2066
2067	ARB_ERROR Aligner::run() {
2068	fa_assert(GB_get_transaction_level(gb_main) == 0); // no open transaction allowed here!
2069	ARB_ERROR error = GB_push_transaction(gb_main);
2070
2071	// if neither 'reference' nor 'get_consensus' are specified -> use next relative for (each) sequence:
2072	bool search_by_pt_server = !reference && !get_consensus;
2073
2074	err_count = 0;
2075	wasNotAllowedToAlign = 0; // incremented for every sequence which has higher protection level (and was not aligned)
2076	species_to_mark.clear();
2077
2078	fa_assert(!reference \|\| !get_consensus); // can't do both modes
2079
2080	if (turnAllowed != FA_TURN_NEVER) {
2081	if ((ali_params.range.is_part()) \|\| !search_by_pt_server) {
2082	// if not selected 'Range/Whole sequence' or not selected 'Reference/Auto search..'
2083	turnAllowed = FA_TURN_NEVER; // then disable mirroring for the current call
2084	}
2085	}
2086
2087	if (!error && !alignment) {
2088	alignment = (GB_CSTR)GBT_get_default_alignment(gb_main); // get default alignment
2089	if (!alignment) error = "No default alignment";
2090	}
2091
2092	if (!error && alignment) {
2093	global_alignmentType = GBT_get_alignment_type(gb_main, alignment);
2094	if (search_by_pt_server) {
2095	GB_alignment_type pt_server_alignmentType = GBT_get_alignment_type(gb_main, relSearch.pt_server_alignment);
2096
2097	if (pt_server_alignmentType != GB_AT_RNA &&
2098	pt_server_alignmentType != GB_AT_DNA) {
2099	error = "pt_servers only support RNA/DNA sequences.\n"
2100	"In the aligner window you may specify a RNA/DNA alignment \n"
2101	"and use a pt_server build on that alignment.";
2102	}
2103	}
2104	}
2105
2106	if (!error) {
2107	GBDATA *gb_species_data = GBT_get_species_data(gb_main);
2108	int max_seq_length = GBT_get_alignment_len(gb_main, alignment);
2109
2110	if (reference) error = alignToExplicitReference(gb_species_data, max_seq_length);
2111	else if (get_consensus) error = alignToConsensus(gb_species_data, max_seq_length);
2112	else error = alignToRelatives(gb_species_data, max_seq_length);
2113	}
2114
2115	ClustalV_exit();
2116	unaligned_bases.clear();
2117
2118	error = GB_end_transaction(gb_main, error); // close global transaction
2119
2120	if (wasNotAllowedToAlign>0) {
2121	const char *mess = GBS_global_string("%i species were not aligned (because of protection level)", wasNotAllowedToAlign);
2122	aw_message(mess);
2123	}
2124
2125	if (err_count) {
2126	aw_message_if(error);
2127	error = GBS_global_string("Aligner produced %i error%c", err_count, err_count==1 ? '\0' : 's');
2128	}
2129
2130	if (error_action != FA_NO_ACTION) {
2131	fa_assert(continue_on_error);
2132
2133	GB_transaction ta(gb_main);
2134	GBT_mark_all(gb_main, 0);
2135	for (GBDATAlist::iterator sp = species_to_mark.begin(); sp != species_to_mark.end(); ++sp) {
2136	GB_write_flag(*sp, 1);
2137	}
2138
2139	const char *whatsMarked = (error_action == FA_MARK_ALIGNED) ? "aligned" : "failed";
2140	size_t markCount = species_to_mark.size();
2141	if (markCount>0) {
2142	const char *msg = GBS_global_string("%zu %s species %s been marked",
2143	markCount,
2144	whatsMarked,
2145	(markCount == 1) ? "has" : "have");
2146	aw_message(msg);
2147	}
2148	}
2149
2150	return error;
2151	}
2152
2153	void FastAligner_start(AW_window aw, const AlignDataAccess data_access) {
2154	AW_root *root = aw->get_root();
2155	char *reference = NULp; // align against next relatives
2156	char *toalign = NULp; // align marked species
2157	ARB_ERROR error = NULp;
2158	int get_consensus = 0;
2159	int pt_server_id = -1;
2160
2161	Aligner_get_first_selected_species get_first_selected_species = NULp;
2162	Aligner_get_next_selected_species get_next_selected_species = NULp;
2163
2164	fa_assert(!island_hopper);
2165	if (root->awar(FA_AWAR_USE_ISLAND_HOPPING)->read_int()) {
2166	island_hopper = new IslandHopping;
2167	if (root->awar(FA_AWAR_USE_SECONDARY)->read_int()) {
2168	char *helix_string = data_access->getHelixString();
2169	if (helix_string) {
2170	island_hopper->set_helix(helix_string);
2171	free(helix_string);
2172	}
2173	else {
2174	error = "Warning: No HELIX found. Can't use secondary structure";
2175	}
2176	}
2177
2178	if (!error) {
2179	island_hopper->set_parameters(root->awar(FA_AWAR_ESTIMATE_BASE_FREQ)->read_int(),
2180	root->awar(FA_AWAR_BASE_FREQ_T)->read_float(),
2181	root->awar(FA_AWAR_BASE_FREQ_C)->read_float(),
2182	root->awar(FA_AWAR_BASE_FREQ_A)->read_float(),
2183	root->awar(FA_AWAR_BASE_FREQ_C)->read_float(),
2184	root->awar(FA_AWAR_SUBST_PARA_CT)->read_float(),
2185	root->awar(FA_AWAR_SUBST_PARA_AT)->read_float(),
2186	root->awar(FA_AWAR_SUBST_PARA_GT)->read_float(),
2187	root->awar(FA_AWAR_SUBST_PARA_AC)->read_float(),
2188	root->awar(FA_AWAR_SUBST_PARA_CG)->read_float(),
2189	root->awar(FA_AWAR_SUBST_PARA_AG)->read_float(),
2190	root->awar(FA_AWAR_EXPECTED_DISTANCE)->read_float(),
2191	root->awar(FA_AWAR_STRUCTURE_SUPPLEMENT)->read_float(),
2192	root->awar(FA_AWAR_GAP_A)->read_float(),
2193	root->awar(FA_AWAR_GAP_B)->read_float(),
2194	root->awar(FA_AWAR_GAP_C)->read_float(),
2195	root->awar(FA_AWAR_THRESHOLD)->read_float()
2196	);
2197	}
2198	}
2199
2200	FA_alignTarget alignWhat = static_cast<FA_alignTarget>(root->awar(FA_AWAR_TO_ALIGN)->read_int());
2201	if (!error) {
2202	switch (alignWhat) {
2203	case FA_CURRENT: { // align current species
2204	toalign = root->awar(AWAR_SPECIES_NAME)->read_string();
2205	break;
2206	}
2207	case FA_MARKED: { // align marked species
2208	break;
2209	}
2210	case FA_SELECTED: { // align selected species
2211	get_first_selected_species = data_access->get_first_selected_species;
2212	get_next_selected_species = data_access->get_next_selected_species;
2213	break;
2214	}
2215	default: {
2216	fa_assert(0);
2217	break;
2218	}
2219	}
2220
2221	switch (static_cast<FA_reference>(root->awar(FA_AWAR_REFERENCE)->read_int())) {
2222	case FA_REF_EXPLICIT: // align against specified species
2223	reference = root->awar(FA_AWAR_REFERENCE_NAME)->read_string();
2224	break;
2225
2226	case FA_REF_CONSENSUS: // align against group consensus
2227	if (data_access->get_group_consensus) {
2228	get_consensus = 1;
2229	}
2230	else {
2231	error = "Can't get group consensus here.";
2232	}
2233	break;
2234
2235	case FA_REF_RELATIVES: // align against species searched via pt_server
2236	pt_server_id = root->awar(AWAR_PT_SERVER)->read_int();
2237	if (pt_server_id<0) {
2238	error = "No pt_server selected";
2239	}
2240	break;
2241
2242	default: fa_assert(0);
2243	break;
2244	}
2245	}
2246
2247	RangeList ranges;
2248	bool autoRestrictRange4nextRelSearch = true;
2249
2250	if (!error) {
2251	switch (static_cast<FA_range>(root->awar(FA_AWAR_RANGE)->read_int())) {
2252	case FA_WHOLE_SEQUENCE:
2253	autoRestrictRange4nextRelSearch = false;
2254	ranges.add(PosRange::whole());
2255	break;
2256
2257	case FA_AROUND_CURSOR: {
2258	int curpos = root->awar(AWAR_CURSOR_POSITION_LOCAL)->read_int();
2259	int size = root->awar(FA_AWAR_AROUND)->read_int();
2260
2261	ranges.add(PosRange(curpos-size, curpos+size));
2262	break;
2263	}
2264	case FA_SELECTED_RANGE: {
2265	PosRange range;
2266	if (!data_access->get_selected_range(range)) {
2267	error = "There is no selected species!";
2268	}
2269	else {
2270	ranges.add(range);
2271	}
2272	break;
2273	}
2274
2275	case FA_SAI_MULTI_RANGE: {
2276	GB_transaction ta(data_access->gb_main);
2277
2278	char *sai_name = root->awar(FA_AWAR_SAI_RANGE_NAME)->read_string();
2279	char *aliuse = GBT_get_default_alignment(data_access->gb_main);
2280
2281	GBDATA *gb_sai = GBT_expect_SAI(data_access->gb_main, sai_name);
2282	if (!gb_sai) error = GB_await_error();
2283	else {
2284	GBDATA *gb_data = GBT_find_sequence(gb_sai, aliuse);
2285	if (!gb_data) {
2286	error = GB_have_error()
2287	? GB_await_error()
2288	: GBS_global_string("SAI '%s' has no data in alignment '%s'", sai_name, aliuse);
2289	}
2290	else {
2291	char *sai_data = GB_read_string(gb_data); // @@@ NOT_ALL_SAI_HAVE_DATA
2292	char *set_bits = root->awar(FA_AWAR_SAI_RANGE_CHARS)->read_string();
2293
2294	ranges = build_RangeList_from_string(sai_data, set_bits, false);
2295
2296	free(set_bits);
2297	free(sai_data);
2298	}
2299	}
2300	free(aliuse);
2301	free(sai_name);
2302	break;
2303	}
2304	}
2305	}
2306
2307	if (!error) {
2308	for (RangeList::iterator r = ranges.begin(); r != ranges.end() && !error; ++r) {
2309	PosRange range = *r;
2310
2311	GBDATA *gb_main = data_access->gb_main;
2312	char *editor_alignment = NULp;
2313	long alignment_length;
2314	{
2315	GB_transaction ta(gb_main);
2316	char *default_alignment = GBT_get_default_alignment(gb_main);
2317
2318	alignment_length = GBT_get_alignment_len(gb_main, default_alignment);
2319
2320	editor_alignment = root->awar_string(AWAR_EDITOR_ALIGNMENT, default_alignment)->read_string();
2321	free(default_alignment);
2322
2323	}
2324
2325	char *pt_server_alignment = root->awar(FA_AWAR_PT_SERVER_ALIGNMENT)->read_string();
2326	PosRange relRange = PosRange::whole(); // unrestricted!
2327
2328	if (autoRestrictRange4nextRelSearch) {
2329	AW_awar *awar_relrange = root->awar(FA_AWAR_RELATIVE_RANGE);
2330	const char *relrange = awar_relrange->read_char_pntr();
2331	if (relrange[0]) {
2332	int region_plus = atoi(relrange);
2333
2334	fa_assert(range.is_part());
2335
2336	relRange = PosRange(range.start()-region_plus, range.end()+region_plus); // restricted
2337	awar_relrange->write_as_string(GBS_global_string("%i", region_plus)); // set awar to detected value (avoid misbehavior when it contains ' ' or similar)
2338	}
2339	}
2340
2341	if (island_hopper) {
2342	island_hopper->set_range(ExplicitRange(range, alignment_length));
2343	range = PosRange::whole();
2344	}
2345
2346	SearchRelativeParams relSearch(new PT_FamilyFinder(gb_main,
2347	pt_server_id,
2348	root->awar(AWAR_NN_OLIGO_LEN)->read_int(),
2349	root->awar(AWAR_NN_MISMATCHES)->read_int(),
2350	root->awar(AWAR_NN_FAST_MODE)->read_int(),
2351	root->awar(AWAR_NN_REL_MATCHES)->read_int(),
2352	RSS_BOTH_MIN), // old scaling as b4 [8520] @@@ make configurable
2353	pt_server_alignment,
2354	root->awar(FA_AWAR_NEXT_RELATIVES)->read_int());
2355
2356	relSearch.getFamilyFinder()->restrict_2_region(relRange);
2357
2358	struct AlignParams ali_params = {
2359	static_cast<FA_report>(root->awar(FA_AWAR_REPORT)->read_int()),
2360	bool(root->awar(FA_AWAR_SHOW_GAPS_MESSAGES)->read_int()),
2361	range
2362	};
2363
2364	{
2365	arb_progress progress("FastAligner");
2366	progress.allow_title_reuse();
2367
2368	int cont_on_error = root->awar(FA_AWAR_CONTINUE_ON_ERROR)->read_int();
2369
2370	Aligner aligner(gb_main,
2371	alignWhat,
2372	editor_alignment,
2373	toalign,
2374	get_first_selected_species,
2375	get_next_selected_species,
2376
2377	reference,
2378	get_consensus ? data_access->get_group_consensus : NULp,
2379	relSearch,
2380
2381	static_cast<FA_turn>(root->awar(FA_AWAR_MIRROR)->read_int()),
2382	ali_params,
2383	root->awar(FA_AWAR_PROTECTION)->read_int(),
2384	cont_on_error,
2385	(FA_errorAction)root->awar(FA_AWAR_ACTION_ON_ERROR)->read_int());
2386	error = aligner.run();
2387
2388	if (error && cont_on_error) {
2389	aw_message_if(error);
2390	error = NULp;
2391	}
2392	}
2393
2394	free(pt_server_alignment);
2395	free(editor_alignment);
2396	}
2397	}
2398
2399	if (island_hopper) {
2400	delete island_hopper;
2401	island_hopper = NULp;
2402	}
2403
2404	if (toalign) free(toalign);
2405	aw_message_if(error);
2406	if (data_access->do_refresh) data_access->refresh_display();
2407	}
2408
2409	static void nullcb() { }
2410
2411	void FastAligner_create_variables(AW_root *root, AW_default db1) {
2412	root->awar_string(FA_AWAR_REFERENCE_NAME, "", db1);
2413
2414	root->awar_int(FA_AWAR_TO_ALIGN, FA_CURRENT, db1);
2415	root->awar_int(FA_AWAR_REFERENCE, FA_REF_EXPLICIT, db1);
2416	root->awar_int(FA_AWAR_RANGE, FA_WHOLE_SEQUENCE, db1);
2417
2418	AW_awar *ali_protect = root->awar_int(FA_AWAR_PROTECTION, 0, db1);
2419	if (ARB_in_novice_mode(root)) {
2420	ali_protect->write_int(0); // reset protection for noobs
2421	}
2422
2423	root->awar_int(FA_AWAR_AROUND, 25, db1);
2424	root->awar_int(FA_AWAR_MIRROR, FA_TURN_INTERACTIVE, db1);
2425	root->awar_int(FA_AWAR_REPORT, FA_NO_REPORT, db1);
2426	root->awar_int(FA_AWAR_SHOW_GAPS_MESSAGES, 1, db1);
2427	root->awar_int(FA_AWAR_CONTINUE_ON_ERROR, 1, db1);
2428	root->awar_int(FA_AWAR_ACTION_ON_ERROR, FA_NO_ACTION, db1);
2429	root->awar_int(FA_AWAR_USE_SECONDARY, 0, db1);
2430	root->awar_int(AWAR_PT_SERVER, 0, db1);
2431	root->awar_int(FA_AWAR_NEXT_RELATIVES, 1, db1)->set_minmax(1, 100);
2432
2433	root->awar_string(FA_AWAR_RELATIVE_RANGE, "", db1);
2434	root->awar_string(FA_AWAR_PT_SERVER_ALIGNMENT, root->awar(AWAR_DEFAULT_ALIGNMENT)->read_char_pntr(), db1);
2435
2436	root->awar_string(FA_AWAR_SAI_RANGE_NAME, "", db1);
2437	root->awar_string(FA_AWAR_SAI_RANGE_CHARS, "x1", db1);
2438
2439	// island hopping:
2440
2441	root->awar_int(FA_AWAR_USE_ISLAND_HOPPING, 0, db1);
2442
2443	root->awar_int(FA_AWAR_ESTIMATE_BASE_FREQ, 1, db1);
2444
2445	root->awar_float(FA_AWAR_BASE_FREQ_A, 0.25, db1);
2446	root->awar_float(FA_AWAR_BASE_FREQ_C, 0.25, db1);
2447	root->awar_float(FA_AWAR_BASE_FREQ_G, 0.25, db1);
2448	root->awar_float(FA_AWAR_BASE_FREQ_T, 0.25, db1);
2449
2450	root->awar_float(FA_AWAR_SUBST_PARA_AC, 1.0, db1);
2451	root->awar_float(FA_AWAR_SUBST_PARA_AG, 4.0, db1);
2452	root->awar_float(FA_AWAR_SUBST_PARA_AT, 1.0, db1);
2453	root->awar_float(FA_AWAR_SUBST_PARA_CG, 1.0, db1);
2454	root->awar_float(FA_AWAR_SUBST_PARA_CT, 4.0, db1);
2455	root->awar_float(FA_AWAR_SUBST_PARA_GT, 1.0, db1);
2456
2457	root->awar_float(FA_AWAR_EXPECTED_DISTANCE, 0.3, db1);
2458	root->awar_float(FA_AWAR_STRUCTURE_SUPPLEMENT, 0.5, db1);
2459	root->awar_float(FA_AWAR_THRESHOLD, 0.005, db1);
2460
2461	root->awar_float(FA_AWAR_GAP_A, 8.0, db1);
2462	root->awar_float(FA_AWAR_GAP_B, 4.0, db1);
2463	root->awar_float(FA_AWAR_GAP_C, 7.0, db1);
2464
2465	AWTC_create_common_next_neighbour_vars(root, makeRootCallback(nullcb));
2466	}
2467
2468	void FastAligner_set_align_current(AW_root *root, AW_default db1) {
2469	root->awar_int(FA_AWAR_TO_ALIGN, FA_CURRENT, db1);
2470	}
2471
2472	void FastAligner_set_reference_species(AW_root *root) {
2473	// sets the aligner reference species to current species
2474	char *specName = root->awar(AWAR_SPECIES_NAME)->read_string();
2475	root->awar(FA_AWAR_REFERENCE_NAME)->write_string(specName);
2476	free(specName);
2477	}
2478
2479	static AW_window create_island_hopping_window(AW_root root) {
2480	AW_window_simple *aws = new AW_window_simple;
2481
2482	aws->init(root, "ISLAND_HOPPING_PARA", "Parameters for Island Hopping");
2483	aws->load_xfig("faligner/islandhopping.fig");
2484
2485	aws->at("close");
2486	aws->callback(AW_POPDOWN);
2487	aws->create_button("CLOSE", "CLOSE", "O");
2488
2489	aws->at("help");
2490	aws->callback(makeHelpCallback("islandhopping.hlp"));
2491	aws->create_button("HELP", "HELP");
2492
2493	aws->at("use_secondary");
2494	aws->label("Use secondary structure (only for re-align)");
2495	aws->create_toggle(FA_AWAR_USE_SECONDARY);
2496
2497	aws->at("freq");
2498	aws->create_toggle_field(FA_AWAR_ESTIMATE_BASE_FREQ, "Base freq.", "B");
2499	aws->insert_default_toggle("Estimate", "E", 1);
2500	aws->insert_toggle("Define here: ", "D", 0);
2501	aws->update_toggle_field();
2502
2503	aws->at("freq_a"); aws->label("A:"); aws->create_input_field(FA_AWAR_BASE_FREQ_A, 4);
2504	aws->at("freq_c"); aws->label("C:"); aws->create_input_field(FA_AWAR_BASE_FREQ_C, 4);
2505	aws->at("freq_g"); aws->label("G:"); aws->create_input_field(FA_AWAR_BASE_FREQ_G, 4);
2506	aws->at("freq_t"); aws->label("T:"); aws->create_input_field(FA_AWAR_BASE_FREQ_T, 4);
2507
2508	int xpos[4], ypos[4];
2509	{
2510	aws->button_length(1);
2511
2512	int dummy;
2513	aws->at("h_a"); aws->get_at_position(&xpos[0], &dummy); aws->create_button(NULp, "A");
2514	aws->at("h_c"); aws->get_at_position(&xpos[1], &dummy); aws->create_button(NULp, "C");
2515	aws->at("h_g"); aws->get_at_position(&xpos[2], &dummy); aws->create_button(NULp, "G");
2516	aws->at("h_t"); aws->get_at_position(&xpos[3], &dummy); aws->create_button(NULp, "T");
2517
2518	aws->at("v_a"); aws->get_at_position(&dummy, &ypos[0]); aws->create_button(NULp, "A");
2519	aws->at("v_c"); aws->get_at_position(&dummy, &ypos[1]); aws->create_button(NULp, "C");
2520	aws->at("v_g"); aws->get_at_position(&dummy, &ypos[2]); aws->create_button(NULp, "G");
2521	aws->at("v_t"); aws->get_at_position(&dummy, &ypos[3]); aws->create_button(NULp, "T");
2522	}
2523
2524	aws->at("subst"); aws->create_button(NULp, "Substitution rate parameters:");
2525
2526	#define XOFF -25
2527	#define YOFF 0
2528
2529	aws->at(xpos[1]+XOFF, ypos[0]+YOFF); aws->create_input_field(FA_AWAR_SUBST_PARA_AC, 4);
2530	aws->at(xpos[2]+XOFF, ypos[0]+YOFF); aws->create_input_field(FA_AWAR_SUBST_PARA_AG, 4);
2531	aws->at(xpos[3]+XOFF, ypos[0]+YOFF); aws->create_input_field(FA_AWAR_SUBST_PARA_AT, 4);
2532	aws->at(xpos[2]+XOFF, ypos[1]+YOFF); aws->create_input_field(FA_AWAR_SUBST_PARA_CG, 4);
2533	aws->at(xpos[3]+XOFF, ypos[1]+YOFF); aws->create_input_field(FA_AWAR_SUBST_PARA_CT, 4);
2534	aws->at(xpos[3]+XOFF, ypos[2]+YOFF); aws->create_input_field(FA_AWAR_SUBST_PARA_GT, 4);
2535
2536	#undef XOFF
2537	#undef YOFF
2538
2539	aws->label_length(22);
2540
2541	aws->at("dist");
2542	aws->label("Expected distance");
2543	aws->create_input_field(FA_AWAR_EXPECTED_DISTANCE, 5);
2544
2545	aws->at("supp");
2546	aws->label("Structure supplement");
2547	aws->create_input_field(FA_AWAR_STRUCTURE_SUPPLEMENT, 5);
2548
2549	aws->at("thres");
2550	aws->label("Threshold");
2551	aws->create_input_field(FA_AWAR_THRESHOLD, 5);
2552
2553	aws->label_length(10);
2554
2555	aws->at("gapA");
2556	aws->label("Gap A");
2557	aws->create_input_field(FA_AWAR_GAP_A, 5);
2558
2559	aws->at("gapB");
2560	aws->label("Gap B");
2561	aws->create_input_field(FA_AWAR_GAP_B, 5);
2562
2563	aws->at("gapC");
2564	aws->label("Gap C");
2565	aws->create_input_field(FA_AWAR_GAP_C, 5);
2566
2567	return aws;
2568	}
2569
2570	static AW_window create_family_settings_window(AW_root root) {
2571	static AW_window_simple *aws = NULp;
2572
2573	if (!aws) {
2574	aws = new AW_window_simple;
2575
2576	aws->init(root, "FAMILY_PARAMS", "Family search parameters");
2577	aws->load_xfig("faligner/family_settings.fig");
2578
2579	aws->at("close");
2580	aws->callback(AW_POPDOWN);
2581	aws->create_button("CLOSE", "CLOSE", "O");
2582
2583	aws->at("help");
2584	aws->callback(makeHelpCallback("next_neighbours_common.hlp"));
2585	aws->create_button("HELP", "HELP");
2586
2587	aws->auto_space(5, 5);
2588	AWTC_create_common_next_neighbour_fields(aws, 200);
2589	}
2590
2591	return aws;
2592	}
2593
2594	static AWT_config_mapping_def aligner_config_mapping[] = {
2595	{ FA_AWAR_USE_ISLAND_HOPPING, "island" },
2596	{ FA_AWAR_TO_ALIGN, "target" },
2597	{ FA_AWAR_REFERENCE, "ref" },
2598	{ FA_AWAR_REFERENCE_NAME, "refname" },
2599	{ FA_AWAR_RELATIVE_RANGE, "relrange" },
2600	{ FA_AWAR_NEXT_RELATIVES, "relatives" },
2601	{ FA_AWAR_PT_SERVER_ALIGNMENT, "ptali" },
2602
2603	// relative-search specific parameters from subwindow (create_family_settings_window)
2604	// same as ../DB_UI/ui_species.cxx@RELATIVES_CONFIG
2605	{ AWAR_NN_OLIGO_LEN, "oligolen" },
2606	{ AWAR_NN_MISMATCHES, "mismatches" },
2607	{ AWAR_NN_FAST_MODE, "fastmode" },
2608	{ AWAR_NN_REL_MATCHES, "relmatches" },
2609	{ AWAR_NN_REL_SCALING, "relscaling" },
2610
2611	// island-hopping parameters (create_island_hopping_window)
2612	{ FA_AWAR_USE_SECONDARY, "use2nd" },
2613	{ FA_AWAR_ESTIMATE_BASE_FREQ, "estbasefreq" },
2614	{ FA_AWAR_BASE_FREQ_A, "freq_a" },
2615	{ FA_AWAR_BASE_FREQ_C, "freq_c" },
2616	{ FA_AWAR_BASE_FREQ_G, "freq_g" },
2617	{ FA_AWAR_BASE_FREQ_T, "freq_t" },
2618	{ FA_AWAR_SUBST_PARA_AC, "subst_ac" },
2619	{ FA_AWAR_SUBST_PARA_AG, "subst_ag" },
2620	{ FA_AWAR_SUBST_PARA_AT, "subst_at" },
2621	{ FA_AWAR_SUBST_PARA_CG, "subst_cg" },
2622	{ FA_AWAR_SUBST_PARA_CT, "subst_ct" },
2623	{ FA_AWAR_SUBST_PARA_GT, "subst_gt" },
2624	{ FA_AWAR_EXPECTED_DISTANCE, "distance" },
2625	{ FA_AWAR_STRUCTURE_SUPPLEMENT, "supplement" },
2626	{ FA_AWAR_THRESHOLD, "threshold" },
2627	{ FA_AWAR_GAP_A, "gap_a" },
2628	{ FA_AWAR_GAP_B, "gap_b" },
2629	{ FA_AWAR_GAP_C, "gap_c" },
2630
2631	{ NULp, NULp }
2632	};
2633
2634	AW_window FastAligner_create_window(AW_root root, const AlignDataAccess *data_access) {
2635	AW_window_simple *aws = new AW_window_simple;
2636
2637	aws->init(root, "INTEGRATED_ALIGNERS", INTEGRATED_ALIGNERS_TITLE);
2638	aws->load_xfig("faligner/faligner.fig");
2639
2640	aws->label_length(10);
2641	aws->button_length(10);
2642
2643	aws->at("close");
2644	aws->callback(AW_POPDOWN);
2645	aws->create_button("CLOSE", "CLOSE", "O");
2646
2647	aws->at("help");
2648	aws->callback(makeHelpCallback("faligner.hlp"));
2649	aws->create_button("HELP", "HELP");
2650
2651	aws->at("aligner");
2652	aws->create_toggle_field(FA_AWAR_USE_ISLAND_HOPPING, "Aligner", "A");
2653	aws->insert_default_toggle("Fast aligner", "F", 0);
2654	aws->sens_mask(AWM_EXP);
2655	aws->insert_toggle ("Island Hopping", "I", 1);
2656	aws->sens_mask(AWM_ALL);
2657	aws->update_toggle_field();
2658
2659	aws->button_length(12);
2660	aws->at("island_para");
2661	aws->callback(create_island_hopping_window);
2662	aws->sens_mask(AWM_EXP);
2663	aws->create_button("island_para", "Parameters", "");
2664	aws->sens_mask(AWM_ALL);
2665
2666	aws->button_length(10);
2667
2668	aws->at("rev_compl");
2669	aws->callback(makeWindowCallback(build_reverse_complement, data_access));
2670	aws->create_button("reverse_complement", "Turn now!", "");
2671
2672	aws->at("what");
2673	aws->create_toggle_field(FA_AWAR_TO_ALIGN, "Align what?", "A");
2674	aws->insert_toggle ("Current Species:", "A", FA_CURRENT);
2675	aws->insert_default_toggle("Marked Species", "M", FA_MARKED);
2676	aws->insert_toggle ("Selected Species", "S", FA_SELECTED);
2677	aws->update_toggle_field();
2678
2679	aws->at("swhat");
2680	aws->create_input_field(AWAR_SPECIES_NAME, 2);
2681
2682	aws->at("against");
2683	aws->create_toggle_field(FA_AWAR_REFERENCE, "Reference", "");
2684	aws->insert_toggle ("Species by name:", "S", FA_REF_EXPLICIT);
2685	aws->insert_toggle ("Group consensus", "K", FA_REF_CONSENSUS);
2686	aws->insert_default_toggle("Auto search by pt_server:", "A", FA_REF_RELATIVES);
2687	aws->update_toggle_field();
2688
2689	aws->at("sagainst");
2690	aws->create_input_field(FA_AWAR_REFERENCE_NAME, 2);
2691
2692	aws->at("copy");
2693	aws->callback(RootAsWindowCallback::simple(FastAligner_set_reference_species));
2694	aws->create_button("Copy", "Copy", "");
2695
2696	aws->label_length(0);
2697	aws->at("pt_server");
2698	awt_create_PTSERVER_selection_button(aws, AWAR_PT_SERVER);
2699
2700	aws->label_length(23);
2701	aws->at("relrange");
2702	aws->label("Data from range only, plus");
2703	aws->create_input_field(FA_AWAR_RELATIVE_RANGE, 3);
2704
2705	aws->at("relatives");
2706	aws->label("Number of relatives to use");
2707	aws->create_input_field(FA_AWAR_NEXT_RELATIVES, 3);
2708
2709	aws->label_length(9);
2710	aws->at("use_ali");
2711	aws->label("Alignment");
2712	aws->create_input_field(FA_AWAR_PT_SERVER_ALIGNMENT, 12);
2713
2714	aws->at("relSett");
2715	aws->callback(create_family_settings_window);
2716	aws->create_autosize_button("Settings", "More settings", "");
2717
2718	// Range
2719
2720	aws->label_length(10);
2721	aws->at("range");
2722	aws->create_toggle_field(FA_AWAR_RANGE, "Range", "");
2723	aws->insert_default_toggle("Whole sequence", "", FA_WHOLE_SEQUENCE);
2724	aws->insert_toggle ("Positions around cursor: ", "", FA_AROUND_CURSOR);
2725	aws->insert_toggle ("Selected range", "", FA_SELECTED_RANGE);
2726	aws->insert_toggle ("Multi-Range by SAI", "", FA_SAI_MULTI_RANGE);
2727	aws->update_toggle_field();
2728
2729	aws->at("around");
2730	aws->create_input_field(FA_AWAR_AROUND, 2);
2731
2732	aws->at("sai");
2733	awt_create_SAI_selection_button(data_access->gb_main, aws, FA_AWAR_SAI_RANGE_NAME);
2734
2735	aws->at("rchars");
2736	aws->create_input_field(FA_AWAR_SAI_RANGE_CHARS, 2);
2737
2738	// Protection
2739
2740	aws->at("protection");
2741	aws->label("Protection");
2742	aws->create_option_menu(FA_AWAR_PROTECTION, true);
2743	aws->insert_default_option("0", NULp, 0);
2744	aws->insert_option ("1", NULp, 1);
2745	aws->insert_option ("2", NULp, 2);
2746	aws->insert_option ("3", NULp, 3);
2747	aws->insert_option ("4", NULp, 4);
2748	aws->insert_option ("5", NULp, 5);
2749	aws->insert_option ("6", NULp, 6);
2750	aws->update_option_menu();
2751
2752	// MirrorCheck
2753
2754	aws->at("mirror");
2755	aws->label("Turn check");
2756	aws->create_option_menu(FA_AWAR_MIRROR, true);
2757	aws->insert_option ("Never turn sequence", "", FA_TURN_NEVER);
2758	aws->insert_default_option("User acknowledgment ", "", FA_TURN_INTERACTIVE);
2759	aws->insert_option ("Automatically turn sequence", "", FA_TURN_ALWAYS);
2760	aws->update_option_menu();
2761
2762	// Report
2763
2764	aws->at("insert");
2765	aws->label("Report");
2766	aws->create_option_menu(FA_AWAR_REPORT, true);
2767	aws->insert_option ("No report", "", FA_NO_REPORT);
2768	aws->sens_mask(AWM_EXP);
2769	aws->insert_default_option("Report to temporary entries", "", FA_TEMP_REPORT);
2770	aws->insert_option ("Report to resident entries", "", FA_REPORT);
2771	aws->sens_mask(AWM_ALL);
2772	aws->update_option_menu();
2773
2774	aws->at("gaps");
2775	aws->create_toggle(FA_AWAR_SHOW_GAPS_MESSAGES);
2776
2777	aws->at("continue");
2778	aws->create_toggle(FA_AWAR_CONTINUE_ON_ERROR);
2779
2780	aws->at("on_failure");
2781	aws->label("On failure");
2782	aws->create_option_menu(FA_AWAR_ACTION_ON_ERROR, true);
2783	aws->insert_default_option("do nothing", "", FA_NO_ACTION);
2784	aws->insert_option ("mark failed", "", FA_MARK_FAILED);
2785	aws->insert_option ("mark aligned", "", FA_MARK_ALIGNED);
2786	aws->update_option_menu();
2787
2788	aws->at("align");
2789	aws->callback(makeWindowCallback(FastAligner_start, data_access));
2790	aws->highlight();
2791	aws->create_button("GO", "GO", "G");
2792
2793	aws->at("config");
2794	AWT_insert_config_manager(aws, AW_ROOT_DEFAULT, "aligner", aligner_config_mapping);
2795
2796	return aws;
2797	}
2798
2799	// --------------------------------------------------------------------------------
2800
2801	#ifdef UNIT_TESTS
2802
2803	#include <test_unit.h>
2804
2805	// ---------------------
2806	// OligoCounter
2807
2808	#include <map>
2809	#include <string>
2810
2811	using std::map;
2812	using std::string;
2813
2814	typedef map<string, size_t> OligoCount;
2815
2816	class OligoCounter {
2817	size_t oligo_len;
2818	size_t datasize;
2819
2820	mutable OligoCount occurrence;
2821
2822	static string removeGaps(const char *seq) {
2823	size_t len = strlen(seq);
2824	string nogaps;
2825	nogaps.reserve(len);
2826
2827	for (size_t p = 0; p<len; ++p) {
2828	char c = seq[p];
2829	if (!is_gap(c)) nogaps.append(1, c);
2830	}
2831	return nogaps;
2832	}
2833
2834	void count_oligos(const string& seq) {
2835	occurrence.clear();
2836	size_t max_pos = seq.length()-oligo_len;
2837	for (size_t p = 0; p <= max_pos; ++p) {
2838	string oligo(seq, p, oligo_len);
2839	occurrence[oligo]++;
2840	}
2841	}
2842
2843	public:
2844	OligoCounter()
2845	: oligo_len(0),
2846	datasize(0)
2847	{}
2848	OligoCounter(const char *seq, size_t oligo_len_)
2849	: oligo_len(oligo_len_)
2850	{
2851	string seq_nogaps = removeGaps(seq);
2852	datasize = seq_nogaps.length();
2853	count_oligos(seq_nogaps);
2854	}
2855
2856	size_t oligo_count(const char *oligo) {
2857	fa_assert(strlen(oligo) == oligo_len);
2858	return occurrence[oligo];
2859	}
2860
2861	size_t similarity_score(const OligoCounter& other) const {
2862	size_t score = 0;
2863	if (oligo_len == other.oligo_len) {
2864	for (OligoCount::const_iterator o = occurrence.begin(); o != occurrence.end(); ++o) {
2865	const string& oligo = o->first;
2866	size_t count = o->second;
2867
2868	score += min(count, other.occurrence[oligo]);
2869	}
2870	}
2871	return score;
2872	}
2873
2874	size_t getDataSize() const { return datasize; }
2875	};
2876
2877	void TEST_OligoCounter() {
2878	OligoCounter oc1("CCAGGT", 3);
2879	OligoCounter oc2("GGTCCA", 3);
2880	OligoCounter oc2_gaps("..GGT--CCA..", 3);
2881	OligoCounter oc3("AGGTCC", 3);
2882	OligoCounter oc4("AGGTCCAGG", 3);
2883
2884	TEST_EXPECT_EQUAL(oc1.oligo_count("CCA"), 1);
2885	TEST_EXPECT_ZERO(oc1.oligo_count("CCG"));
2886	TEST_EXPECT_EQUAL(oc4.oligo_count("AGG"), 2);
2887
2888	int sc1_2 = oc1.similarity_score(oc2);
2889	int sc2_1 = oc2.similarity_score(oc1);
2890	TEST_EXPECT_EQUAL(sc1_2, sc2_1);
2891
2892	int sc1_2gaps = oc1.similarity_score(oc2_gaps);
2893	TEST_EXPECT_EQUAL(sc1_2, sc1_2gaps);
2894
2895	int sc1_3 = oc1.similarity_score(oc3);
2896	int sc2_3 = oc2.similarity_score(oc3);
2897	int sc3_4 = oc3.similarity_score(oc4);
2898
2899	TEST_EXPECT_EQUAL(sc1_2, 2); // common oligos (CCA GGT)
2900	TEST_EXPECT_EQUAL(sc1_3, 2); // common oligos (AGG GGT)
2901	TEST_EXPECT_EQUAL(sc2_3, 3); // common oligos (GGT GTC TCC)
2902
2903	TEST_EXPECT_EQUAL(sc3_4, 4);
2904	}
2905
2906	// -------------------------
2907	// FakeFamilyFinder
2908
2909	class FakeFamilyFinder: public FamilyFinder { // derived from a Noncopyable
2910	// used by unit tests to detect next relatives instead of asking the pt-server
2911
2912	GBDATA *gb_main;
2913	string ali_name;
2914	map<string, OligoCounter> oligos_counted; // key = species name
2915	PosRange counted_for_range;
2916	size_t oligo_len;
2917
2918	public:
2919	FakeFamilyFinder(GBDATA *gb_main_, string ali_name_, bool rel_matches_, size_t oligo_len_)
2920	: FamilyFinder(rel_matches_, RSS_BOTH_MIN),
2921	gb_main(gb_main_),
2922	ali_name(ali_name_),
2923	counted_for_range(PosRange::whole()),
2924	oligo_len(oligo_len_)
2925	{}
2926
2927	GB_ERROR searchFamily(const char *sequence, FF_complement compl_mode, int max_results, double min_score) OVERRIDE {
2928	// 'sequence' has to contain full sequence or part corresponding to 'range'
2929
2930	TEST_EXPECT_EQUAL(compl_mode, FF_FORWARD); // not fit for other modes
2931
2932	delete_family_list();
2933
2934	OligoCounter seq_oligo_count(sequence, oligo_len);
2935
2936	if (range != counted_for_range) {
2937	oligos_counted.clear(); // forget results for different range
2938	counted_for_range = range;
2939	}
2940
2941	char *buffer = NULp;
2942	int buffersize = 0;
2943
2944	bool partial_match = range.is_part();
2945
2946	GB_transaction ta(gb_main);
2947	int results = 0;
2948
2949	for (GBDATA *gb_species = GBT_first_species(gb_main);
2950	gb_species && results<max_results;
2951	gb_species = GBT_next_species(gb_species))
2952	{
2953	string name = GBT_get_name_or_description(gb_species);
2954	if (oligos_counted.find(name) == oligos_counted.end()) {
2955	GBDATA *gb_data = GBT_find_sequence(gb_species, ali_name.c_str());
2956	const char *spec_seq = GB_read_char_pntr(gb_data);
2957
2958	if (partial_match) {
2959	int spec_seq_len = GB_read_count(gb_data);
2960	int range_len = ExplicitRange(range, spec_seq_len).size();
2961
2962	if (buffersize<range_len) {
2963	delete [] buffer;
2964	buffersize = range_len;
2965	buffer = new char[buffersize+1];
2966	}
2967
2968	range.copy_corresponding_part(buffer, spec_seq, spec_seq_len);
2969	oligos_counted[name] = OligoCounter(buffer, oligo_len);
2970	}
2971	else {
2972	oligos_counted[name] = OligoCounter(spec_seq, oligo_len);
2973	}
2974	}
2975
2976	const OligoCounter& spec_oligo_count = oligos_counted[name];
2977	size_t score = seq_oligo_count.similarity_score(spec_oligo_count);
2978
2979	if (score>=min_score) {
2980	FamilyList *newMember = new FamilyList;
2981
2982	newMember->name = strdup(name.c_str());
2983	newMember->matches = score;
2984	newMember->rel_matches = score/spec_oligo_count.getDataSize();
2985	newMember->next = NULp;
2986
2987	family_list = newMember->insertSortedBy_matches(family_list);
2988	results++;
2989	}
2990	}
2991
2992	delete [] buffer;
2993
2994	return NULp;
2995	}
2996	};
2997
2998	// ----------------------------
2999	// test_alignment_data
3000
3001	#include <arb_unit_test.h>
3002
3003	static const char *test_aliname = "ali_test";
3004
3005	static const char get_aligned_data_of(GBDATA gb_main, const char *species_name) {
3006	GB_transaction ta(gb_main);
3007	ARB_ERROR error;
3008	const char *data = NULp;
3009
3010	GBDATA *gb_species = GBT_find_species(gb_main, species_name);
3011	if (!gb_species) error = GB_await_error();
3012	else {
3013	GBDATA *gb_data = GBT_find_sequence(gb_species, test_aliname);
3014	if (!gb_data) error = GB_await_error();
3015	else {
3016	data = GB_read_char_pntr(gb_data);
3017	if (!data) error = GB_await_error();
3018	}
3019	}
3020
3021	TEST_EXPECT_NULL(error.deliver());
3022
3023	return data;
3024	}
3025
3026	static const char get_used_rels_for(GBDATA gb_main, const char *species_name) {
3027	GB_transaction ta(gb_main);
3028	const char *result = NULp;
3029	GBDATA *gb_species = GBT_find_species(gb_main, species_name);
3030	if (!gb_species) result = GBS_global_string("<No such species '%s'>", species_name);
3031	else {
3032	GBDATA *gb_used_rels = GB_search(gb_species, "used_rels", GB_FIND);
3033	if (!gb_used_rels) result = "<No such field 'used_rels'>";
3034	else result = GB_read_char_pntr(gb_used_rels);
3035	}
3036	return result;
3037	}
3038
3039	static GB_ERROR forget_used_relatives(GBDATA *gb_main) {
3040	GB_ERROR error = NULp;
3041	GB_transaction ta(gb_main);
3042	for (GBDATA *gb_species = GBT_first_species(gb_main);
3043	gb_species && !error;
3044	gb_species = GBT_next_species(gb_species))
3045	{
3046	GBDATA *gb_used_rels = GB_search(gb_species, "used_rels", GB_FIND);
3047	if (gb_used_rels) {
3048	error = GB_delete(gb_used_rels);
3049	}
3050	}
3051	return error;
3052	}
3053
3054
3055	#define ALIGNED_DATA_OF(name) get_aligned_data_of(gb_main, name)
3056	#define USED_RELS_FOR(name) get_used_rels_for(gb_main, name)
3057
3058	// ----------------------------------------
3059
3060	static GBDATA *selection_fake_gb_main = NULp;
3061	static GBDATA *selection_fake_gb_last = NULp;
3062
3063	static GBDATA fake_first_selected(int count) {
3064	selection_fake_gb_last = NULp;
3065	*count = 2; // we fake two species as selected ("c1" and "c2")
3066	return GBT_find_species(selection_fake_gb_main, "c1");
3067	}
3068	static GBDATA *fake_next_selected() {
3069	if (!selection_fake_gb_last) {
3070	selection_fake_gb_last = GBT_find_species(selection_fake_gb_main, "c2");
3071	}
3072	else {
3073	selection_fake_gb_last = NULp;
3074	}
3075	return selection_fake_gb_last;
3076	}
3077
3078	static char fake_get_consensus(const char, PosRange range) {
3079	const char *data = get_aligned_data_of(selection_fake_gb_main, "s1");
3080	if (range.is_whole()) return strdup(data);
3081	return ARB_strpartdup(data+range.start(), data+range.end());
3082	}
3083
3084	static void test_install_fakes(GBDATA *gb_main) {
3085	selection_fake_gb_main = gb_main;
3086	}
3087
3088
3089	// ----------------------------------------
3090
3091	static AlignParams test_ali_params = {
3092	FA_NO_REPORT,
3093	false, // showGapsMessages
3094	PosRange::whole()
3095	};
3096
3097	static AlignParams test_ali_params_partial = {
3098	FA_NO_REPORT,
3099	false, // showGapsMessages
3100	PosRange(25, 60)
3101	};
3102
3103	// ----------------------------------------
3104
3105	static struct arb_unit_test::test_alignment_data TestAlignmentData_TargetAndReferenceHandling[] = {
3106	{ 0, "s1", ".........A--UCU-C------C-U-AAACC-CA-A-C-C-G-UAG-UUC--------GAA-U-UGAGG-AC--U-GUAA-CU-C..........." }, // reference
3107	{ 0, "s2", "AUCUCCUAAACCCAACCGUAGUUCGAAUUGAGGACUGUAACUC......................................................" }, // align single sequence (same data as reference)
3108	{ 1, "m1", "UAGAGGAUUUGGGUUGGCAUCAAGCUUAACUCCUGACAUUGAG......................................................" }, // align marked sequences.. (complement of reference)
3109	{ 1, "m2", "...UCCUAAACCAACCCGUAGUUCGAAUUGAGGACUGUAA........................................................." },
3110	{ 1, "m3", "AUC---UAAACCAACCCGUAGUUCGAAUUGAGGACUG---CUC......................................................" },
3111	{ 0, "c1", "AUCUCCUAAACCCAACC--------AAUUGAGGACUGUAACUC......................................................" }, // align selected sequences..
3112	{ 0, "c2", "AUCUCCU------AACCGUAGUUCCCCGAA------ACUGUAACUC..................................................." },
3113	{ 0, "r1", "GAGUUACAGUCCUCAAUUCGGGGAACUACGGUUGGGUUUAGGAGAU..................................................." }, // align by faked pt_server
3114	};
3115
3116	void TEST_Aligner_TargetAndReferenceHandling() {
3117	// performs some alignments to check logic of target and reference handling
3118
3119	GB_shell shell;
3120	ARB_ERROR error;
3121	GBDATA *gb_main = TEST_CREATE_DB(error, test_aliname, TestAlignmentData_TargetAndReferenceHandling, false);
3122
3123	TEST_EXPECT_NULL(error.deliver());
3124
3125	SearchRelativeParams search_relative_params(new FakeFamilyFinder(gb_main, test_aliname, false, 8),
3126	test_aliname,
3127	2);
3128
3129	test_install_fakes(gb_main);
3130	arb_suppress_progress silence;
3131
3132	// bool cont_on_err = true;
3133	bool cont_on_err = false;
3134
3135	TEST_EXPECT_EQUAL(GBT_count_marked_species(gb_main), 3); // we got 3 marked species
3136	{
3137	Aligner aligner(gb_main,
3138	FA_CURRENT,
3139	test_aliname,
3140	"s2", // toalign
3141	NULp, // get_first_selected_species
3142	NULp, // get_next_selected_species
3143	"s1", // reference species
3144	NULp, // get_consensus
3145	search_relative_params, // relative search
3146	FA_TURN_ALWAYS,
3147	test_ali_params,
3148	0,
3149	cont_on_err,
3150	FA_NO_ACTION);
3151	error = aligner.run();
3152	TEST_EXPECT_NULL(error.deliver());
3153	}
3154	TEST_EXPECT_EQUAL(GBT_count_marked_species(gb_main), 3); // we still got 3 marked species
3155	{
3156	Aligner aligner(gb_main,
3157	FA_MARKED,
3158	test_aliname,
3159	NULp, // toalign
3160	NULp, // get_first_selected_species
3161	NULp, // get_next_selected_species
3162	"s1", // reference species
3163	NULp, // get_consensus
3164	search_relative_params, // relative search
3165	FA_TURN_ALWAYS,
3166	test_ali_params,
3167	0,
3168	cont_on_err,
3169	FA_MARK_FAILED);
3170	error = aligner.run();
3171	TEST_EXPECT_NULL(error.deliver());
3172
3173	TEST_EXPECT(!cont_on_err \|\| GBT_count_marked_species(gb_main) == 0); // FA_MARK_FAILED (none failed -> none marked)
3174	}
3175	{
3176	Aligner aligner(gb_main,
3177	FA_SELECTED,
3178	test_aliname,
3179	NULp, // toalign
3180	fake_first_selected, // get_first_selected_species
3181	fake_next_selected, // get_next_selected_species
3182	NULp, // reference species
3183	fake_get_consensus, // get_consensus
3184	search_relative_params, // relative search
3185	FA_TURN_ALWAYS,
3186	test_ali_params,
3187	0,
3188	cont_on_err,
3189	FA_MARK_ALIGNED);
3190	error = aligner.run();
3191	TEST_EXPECT_NULL(error.deliver());
3192
3193	TEST_EXPECT(!cont_on_err \|\| GBT_count_marked_species(gb_main) == 2); // FA_MARK_ALIGNED (2 selected were aligned)
3194	}
3195	{
3196	Aligner aligner(gb_main,
3197	FA_CURRENT,
3198	test_aliname,
3199	"r1", // toalign
3200	NULp, // get_first_selected_species
3201	NULp, // get_next_selected_species
3202	NULp, // reference species
3203	NULp, // get_consensus
3204	search_relative_params, // relative search
3205	FA_TURN_ALWAYS,
3206	test_ali_params,
3207	0,
3208	cont_on_err,
3209	FA_MARK_ALIGNED);
3210
3211	error = aligner.run();
3212	TEST_EXPECT_NULL(error.deliver());
3213
3214	TEST_EXPECT(!cont_on_err \|\| GBT_count_marked_species(gb_main) == 1);
3215	}
3216
3217	TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("s2"), ".........A--UCU-C------C-U-AAACC-CA-A-C-C-G-UAG-UUC--------GAA-U-UGAGG-AC--U-GUAA-CU-C...........");
3218
3219	TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("m1"), ".......UAG--AGG-A------U-U-UGGGU-UG-G-C-A-U-CAA-GCU--------UAA-C-UCCUG-AC--A-UUGAG...............");
3220	TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("m2"), "..............U-C------C-U-AAACC-AA-C-C-C-G-UAG-UUC--------GAA-U-UGAGG-AC--U-GUAA................");
3221	TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("m3"), ".........A--U----------C-U-AAACC-AA-C-C-C-G-UAG-UUC--------GAA-U-UGAGG-AC--U-G----CU-C...........");
3222
3223	TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("c1"), ".........A--UCU-C------C-U-AAACC-CA-A-C-C-------------------AA-U-UGAGG-AC--U-GUAA-CU-C...........");
3224	TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("c2"), ".........A--UCU-C------C-U-AA---------C-C-G-UAG-UUC------------C-CCGAA-AC--U-GUAA-CU-C...........");
3225
3226	TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("r1"), ".........A--UCU-C------C-U-AAACC-CA-A-C-C-G-UAG-UUCCCC-----GAA-U-UGAGG-AC--U-GUAA-CU-C..........."); // here sequence shall be turned!
3227
3228
3229	TEST_EXPECT_EQUAL(USED_RELS_FOR("r1"), "s2:43, s1:1");
3230
3231	// ----------------------------------------------
3232	// now align all others vs next relative
3233
3234	search_relative_params.maxRelatives = 5;
3235	TEST_EXPECT_NO_ERROR(forget_used_relatives(gb_main));
3236
3237	int species_count = ARRAY_ELEMS(TestAlignmentData_TargetAndReferenceHandling);
3238	for (int sp = 0; sp<species_count; ++sp) {
3239	const char *name = TestAlignmentData_TargetAndReferenceHandling[sp].name;
3240	if (strcmp(name, "r1") != 0) { // skip "r1" (already done above)
3241	Aligner aligner(gb_main,
3242	FA_CURRENT,
3243	test_aliname,
3244	name, // toalign
3245	NULp, // get_first_selected_species
3246	NULp, // get_next_selected_species
3247	NULp, // reference species
3248	NULp, // get_consensus
3249	search_relative_params, // relative search
3250	FA_TURN_ALWAYS,
3251	test_ali_params,
3252	0,
3253	cont_on_err,
3254	FA_MARK_ALIGNED);
3255
3256	error = aligner.run();
3257	TEST_EXPECT_NULL(error.deliver());
3258
3259	TEST_EXPECT(!cont_on_err \|\| GBT_count_marked_species(gb_main) == 1);
3260	}
3261	}
3262
3263	TEST_EXPECT_EQUAL(USED_RELS_FOR("s1"), "s2");
3264	TEST_EXPECT_EQUAL(USED_RELS_FOR("s2"), "s1"); // same as done manually
3265	TEST_EXPECT_EQUAL(USED_RELS_FOR("m1"), "r1:42");
3266	TEST_EXPECT_EQUAL(USED_RELS_FOR("m2"), "m3");
3267	TEST_EXPECT_EQUAL(USED_RELS_FOR("m3"), "m2");
3268	TEST_EXPECT_EQUAL(USED_RELS_FOR("c1"), "r1");
3269	TEST_EXPECT_EQUAL(USED_RELS_FOR("c2"), "r1");
3270
3271	// range aligned below (see test_ali_params_partial)
3272	// "-------------------------XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX------------------------------------"
3273	TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("s1"), ".........A--UCU-C------C-U-AAACC-CA-A-C-C-G-UAG-UUC--------GAA-U-UGAGG-AC--U-GUAA-CU-C..........."); // 1st aligning of 's1'
3274	TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("s2"), ".........A--UCU-C------C-U-AAACC-CA-A-C-C-G-UAG-UUC--------GAA-U-UGAGG-AC--U-GUAA-CU-C..........."); // same_as_above (again aligned vs 's1')
3275
3276	TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("m1"), ".........U--AGA-G------G---AUUUG-GG-U-U-G-G-CAU-CAAGCU-----UAA-C-UCCUG-AC--A-UUGAG---------------"); // changed; @@@ bug: no dots at end
3277
3278	TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("m2"), ".........U--C----------C-U-AAACC-AA-C-C-C-G-UAG-UUC--------GAA-U-UGAGG-AC--U-G----UA-A..........."); // changed (1st align vs 's1', this align vs 'm3')
3279	TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("m3"), ".........A--U----------C-U-AAACC-AA-C-C-C-G-UAG-UUC--------GAA-U-UGAGG-AC--U-G----CU-C..........."); // same_as_above (1st align vs 's1', this align vs 'm2')
3280	TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("c1"), ".........A--UCU-C------C-U-AAACC-CA-A-C-C-------------------AA-U-UGAGG-AC--U-GUAA-CU-C..........."); // same_as_above
3281	TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("c2"), ".........A--UCU-C------C-U--------A-A-C-C-G-UAG-UUCCCC-----GA--------A-AC--U-GUAA-CU-C..........."); // changed
3282
3283	// --------------------------------------
3284	// test partial relative search
3285
3286	search_relative_params.getFamilyFinder()->restrict_2_region(test_ali_params_partial.range);
3287	TEST_EXPECT_NO_ERROR(forget_used_relatives(gb_main));
3288
3289	for (int sp = 0; sp<species_count; ++sp) {
3290	const char *name = TestAlignmentData_TargetAndReferenceHandling[sp].name;
3291	Aligner aligner(gb_main,
3292	FA_CURRENT,
3293	test_aliname,
3294	name, // toalign
3295	NULp, // get_first_selected_species
3296	NULp, // get_next_selected_species
3297	NULp, // reference species
3298	NULp, // get_consensus
3299	search_relative_params, // relative search
3300	FA_TURN_NEVER,
3301	test_ali_params_partial,
3302	0,
3303	cont_on_err,
3304	FA_MARK_ALIGNED);
3305
3306	error = aligner.run();
3307	TEST_EXPECT_NULL(error.deliver());
3308
3309	TEST_EXPECT(!cont_on_err \|\| GBT_count_marked_species(gb_main) == 1);
3310	}
3311
3312	TEST_EXPECT_EQUAL(USED_RELS_FOR("s1"), "s2");
3313	TEST_EXPECT_EQUAL(USED_RELS_FOR("s2"), "s1");
3314	TEST_EXPECT_EQUAL(USED_RELS_FOR("m1"), "r1"); // (not really differs)
3315	TEST_EXPECT_EQUAL(USED_RELS_FOR("m2"), "m3");
3316	TEST_EXPECT_EQUAL(USED_RELS_FOR("m3"), "m2");
3317	TEST_EXPECT_EQUAL(USED_RELS_FOR("c1"), "r1");
3318	TEST_EXPECT_EQUAL(USED_RELS_FOR("c2"), "r1");
3319	TEST_EXPECT_EQUAL(USED_RELS_FOR("r1"), "s2:20, c2:3");
3320
3321	// aligned range (see test_ali_params_partial)
3322	// "-------------------------XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX------------------------------------"
3323	TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("s1"), ".........A--UCU-C------C-U-AAACC-CA-A-C-C-G-UAG-UUC--------GAA-U-UGAGG-AC--U-GUAA-CU-C..........."); // same_as_above
3324	TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("s2"), ".........A--UCU-C------C-U-AAACC-CA-A-C-C-G-UAG-UUC--------GAA-U-UGAGG-AC--U-GUAA-CU-C..........."); // same_as_above
3325
3326	TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("m1"), ".........U--AGA-G------G-A-UU-UG-GG-U-U-G-G-CAU-CAAGCU-----UAA-C-UCCUG-AC--A-UUGAG---------------"); // changed
3327	TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("m2"), ".........U--C----------C-U-AAACC-AA-C-C-C-G-UAG-UUC--------GAA-U-UGAGG-AC--U-G----UA-A..........."); // same_as_above
3328	TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("m3"), ".........A--U----------C-U-AAACC-AA-C-C-C-G-UAG-UUC--------GAA-U-UGAGG-AC--U-G----CU-C..........."); // same_as_above
3329
3330	TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("c1"), ".........A--UCU-C------C-U-AAACC-CA-A-C-C-------------------AA-U-UGAGG-AC--U-GUAA-CU-C..........."); // same_as_above
3331	TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("c2"), ".........A--UCU-C------C---------UA-A-C-C-G-UAG-UUCCCC-----GA--------A-AC--U-GUAA-CU-C..........."); // changed
3332
3333	TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("r1"), ".........A--UCU-C------C-U-AAACC-CA-A-C-C-G-UAG-UUCCCC-----GAA-U-UGAGG-AC--U-GUAA-CU-C..........."); // same_as_above
3334
3335	GB_close(gb_main);
3336	}
3337
3338	// ----------------------------------------
3339
3340	static struct arb_unit_test::test_alignment_data TestAlignmentData_checksumError[] = {
3341	{ 0, "MtnK1722", "...G-GGC-C-G............CCC-GG--------CAAUGGGGGCGGCCCGGCGGAC----GG--C-UCAGU-A---AAG-UCGUAACAA-GG-UAG-CCGU-AGGGGAA-CCUG-CGGC-UGGAUCACCUCC....." }, // gets aligned
3342	{ 0, "MhnFormi", "...A-CGA-U-C------------CUUCGG--------GGUCG-U-GG-C-GU-A--C------GG--C-UCAGU-A---AAG-UCGUAACAA-GG-UAG-CCGU-AGGGGAA-CCUG-CGGC-UGGAUCACCUCCU...." }, // 1st relative
3343	{ 0, "MhnT1916", "...A-CGA-A-C------------CUU-GU--------GUUCG-U-GG-C-GA-A--C------GG--C-UCAGU-A---AAG-UCGUAACAA-GG-UAG-CCGU-AGGGGAA-CCUG-CGGC-UGGAUCACCUCCU...." }, // next relative
3344	{ 0, "MthVanni", "...U-GGU-U-U------------C-------------GGCCA-U-GG-C-GG-A--C------GG--C-UCAUU-A---AAG-UCGUAACAA-GG-UAG-CCGU-AGGGGAA-CCUG-CGGC-UGGAUCACCUCC....." }, // next relative
3345	{ 0, "ThcCeler", "...G-GGG-C-G...CC-U---U--------GC--G--CGCAC-C-GG-C-GG-A--C------GG--C-UCAGU-A---AAG-UCGUAACAA-GG-UAG-CCGU-AGGGGAA-CCUA-CGGC-UCGAUCACCUCCU...." }, // next relative
3346	};
3347
3348	void TEST_SLOW_Aligner_checksumError() {
3349	// @@@ SLOW cause this often gets terminated in nightly builds
3350	// no idea why. normally it runs a few ms
3351
3352	// this produces an internal aligner error
3353
3354	GB_shell shell;
3355	ARB_ERROR error;
3356	GBDATA *gb_main = TEST_CREATE_DB(error, test_aliname, TestAlignmentData_checksumError, false);
3357
3358	SearchRelativeParams search_relative_params(new FakeFamilyFinder(gb_main, test_aliname, false, 8),
3359	test_aliname,
3360	10); // use up to 10 relatives
3361
3362	test_install_fakes(gb_main);
3363	arb_suppress_progress silence;
3364
3365	bool cont_on_err = true;
3366	if (!error) {
3367	Aligner aligner(gb_main,
3368	FA_CURRENT,
3369	test_aliname,
3370	"MtnK1722", // toalign
3371	NULp, // get_first_selected_species
3372	NULp, // get_next_selected_species
3373	NULp, // reference species
3374	NULp, // get_consensus
3375	search_relative_params, // relative search
3376	FA_TURN_ALWAYS,
3377	test_ali_params,
3378	0,
3379	cont_on_err,
3380	FA_MARK_ALIGNED);
3381
3382	error = aligner.run();
3383	}
3384	{
3385	GB_ERROR err = error.deliver();
3386	TEST_EXPECT_NULL__BROKEN(err, "Aligner produced 1 error");
3387	}
3388	TEST_EXPECT_EQUAL__BROKEN(USED_RELS_FOR("MtnK1722"), "???", "<No such field 'used_rels'>"); // subsequent failure
3389
3390	GB_close(gb_main);
3391	}
3392
3393	static const char *asstr(LooseBases& ub) {
3394	LooseBases tmp;
3395	while (!ub.is_empty()) tmp.memorize(ub.recall());
3396
3397	const char *result = "";
3398	while (!tmp.is_empty()) {
3399	ExplicitRange r = tmp.recall();
3400	result = GBS_global_string("%s %i/%i", result, r.start(), r.end());
3401	ub.memorize(r);
3402	}
3403	return result;
3404	}
3405
3406	void TEST_BASIC_UnalignedBases() {
3407	LooseBases ub;
3408	TEST_EXPECT(ub.is_empty());
3409	TEST_EXPECT_EQUAL(asstr(ub), "");
3410
3411	// test add+remove
3412	ub.memorize(ExplicitRange(5, 7));
3413	TEST_REJECT(ub.is_empty());
3414	TEST_EXPECT_EQUAL(asstr(ub), " 5/7");
3415
3416	TEST_EXPECT(ub.recall() == ExplicitRange(5, 7));
3417	TEST_EXPECT(ub.is_empty());
3418
3419	ub.memorize(ExplicitRange(2, 4));
3420	TEST_EXPECT_EQUAL(asstr(ub), " 2/4");
3421
3422	ub.memorize(ExplicitRange(4, 9));
3423	TEST_EXPECT_EQUAL(asstr(ub), " 2/4 4/9");
3424
3425	ub.memorize(ExplicitRange(8, 10));
3426	ub.memorize(ExplicitRange(11, 14));
3427	ub.memorize(ExplicitRange(12, 17));
3428	TEST_EXPECT_EQUAL(asstr(ub), " 2/4 4/9 8/10 11/14 12/17");
3429	TEST_EXPECT_EQUAL(asstr(ub), " 2/4 4/9 8/10 11/14 12/17"); // check asstr has no side-effect
3430
3431	{
3432	LooseBases toaddNrecalc;
3433
3434	CompactedSubSequence Old("ACGTACGT", 8, "name1");
3435	CompactedSubSequence New("--A-C--G-T--A-C-G-T", 19, "name2");
3436	// ---------------------- 0123456789012345678
3437
3438	toaddNrecalc.memorize(ExplicitRange(1, 7));
3439	toaddNrecalc.memorize(ExplicitRange(3, 5));
3440	TEST_EXPECT_EQUAL(asstr(toaddNrecalc), " 1/7 3/5");
3441
3442	ub.follow_ali_change_and_append(toaddNrecalc, AliChange(Old, New));
3443
3444	TEST_EXPECT_EQUAL(asstr(ub), " 3/18 8/15 2/4 4/9 8/10 11/14 12/17");
3445	TEST_EXPECT(toaddNrecalc.is_empty());
3446
3447	LooseBases selfRecalc;
3448	selfRecalc.follow_ali_change_and_append(ub, AliChange(New, New));
3449	TEST_EXPECT_EQUAL__BROKEN(asstr(selfRecalc),
3450	" 3/18 8/15 0/6 3/11 8/11 10/15 10/17", // wanted behavior?
3451	" 3/18 8/17 0/6 3/11 8/13 10/15 10/18"); // doc wrong behavior @@@ "8/17", "8/13", "10/18" are wrong
3452
3453	ub.follow_ali_change_and_append(selfRecalc, AliChange(New, Old));
3454	TEST_EXPECT_EQUAL__BROKEN(asstr(ub),
3455	" 1/7 3/5 0/1 1/3 3/3 4/5 4/6", // wanted behavior? (from wanted behavior above)
3456	" 1/7 3/7 0/2 1/4 3/5 4/6 4/7"); // document wrong behavior
3457	TEST_EXPECT_EQUAL__BROKEN(asstr(ub),
3458	" 1/7 3/6 0/1 1/3 3/4 4/5 4/7", // wanted behavior? (from wrong result above)
3459	" 1/7 3/7 0/2 1/4 3/5 4/6 4/7"); // document wrong behavior
3460	}
3461	}
3462
3463
3464	#endif // UNIT_TESTS
3465
3466

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

Download in other formats: