Context Navigation

source: tags/ms_r16q2/SL/FAST_ALIGNER/fast_aligner.cxx

Visit:

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

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

Download in other formats: