Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

fast_aligner.cxx

Visit:

Last change on this file was 18125, checked in by westram, 5 years ago
full update from 'trunk' into 'nameserver' adds: log:branches/cleanup@17627:17650,18105:18112 log:branches/clone@17813:17835 log:branches/config@17733:17737 log:branches/dbui@17881:17889 log:branches/fix@17597:17617 log:branches/fts@17653:18088,18090:18094 log:branches/gcc@17850:17876 log:branches/gui@17754:18092 log:branches/mix@17631:17633,17637:17640 log:trunk@17594:18123
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 136.4 KB

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

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

Context Navigation

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

Download in other formats: