Context Navigation

source: branches/help/DIST/DI_matr.cxx

Visit:

Last change on this file was 18634, checked in by westram, 4 years ago
reintegrates 'root' into 'trunk' implements #449 adds: log:branches/root@18451:18633
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 65.8 KB

Line
1	// =============================================================== //
2	// //
3	// File : DI_matr.cxx //
4	// Purpose : //
5	// //
6	// Institute of Microbiology (Technical University Munich) //
7	// http://www.arb-home.de/ //
8	// //
9	// =============================================================== //
10
11	#include "di_protdist.hxx"
12	#include "di_clusters.hxx"
13	#include "di_view_matrix.hxx"
14	#include "di_awars.hxx"
15
16	#include <neighbourjoin.hxx>
17	#include <AP_filter.hxx>
18	#include <CT_ctree.hxx>
19	#include <ColumnStat.hxx>
20
21	#include <awt.hxx>
22	#include <awt_sel_boxes.hxx>
23	#include <awt_filter.hxx>
24
25	#include <aw_preset.hxx>
26	#include <aw_awars.hxx>
27	#include <aw_file.hxx>
28	#include <aw_msg.hxx>
29	#include <aw_root.hxx>
30
31	#include <gui_aliview.hxx>
32
33	#include <climits>
34	#include <ctime>
35	#include <cmath>
36	#include <arb_sort.h>
37	#include <arb_global_defs.h>
38	#include <macros.hxx>
39	#include <ad_cb.h>
40	#include <awt_TreeAwars.hxx>
41	#include <arb_defs.h>
42
43	using std::string;
44
45	// --------------------------------------------------------------------------------
46
47	#define AWAR_DIST_BOOTSTRAP_COUNT AWAR_DIST_PREFIX "bootstrap/count"
48	#define AWAR_DIST_CANCEL_CHARS AWAR_DIST_PREFIX "cancel/chars"
49	#define AWAR_DIST_MATRIX_AUTO_RECALC AWAR_DIST_PREFIX "recalc"
50
51	#define AWAR_DIST_COLUMN_STAT_NAME AWAR_DIST_COLUMN_STAT_PREFIX "name"
52
53	#define AWAR_DIST_TREE_SORT_NAME "tmp/" AWAR_DIST_TREE_PREFIX "sort_tree_name"
54	#define AWAR_DIST_TREE_COMP_NAME "tmp/" AWAR_DIST_TREE_PREFIX "compr_tree_name"
55	#define AWAR_DIST_TREE_STD_NAME AWAR_DIST_TREE_PREFIX "tree_name"
56	#define AWAR_DIST_MATRIX_AUTO_CALC_TREE AWAR_DIST_TREE_PREFIX "autocalc"
57
58	#define AWAR_DIST_SAVE_MATRIX_TYPE AWAR_DIST_SAVE_MATRIX_BASE "/type"
59	#define AWAR_DIST_SAVE_MATRIX_FILENAME AWAR_DIST_SAVE_MATRIX_BASE "/file_name"
60
61	// --------------------------------------------------------------------------------
62
63	DI_GLOBAL_MATRIX GLOBAL_MATRIX;
64
65	static MatrixDisplay *matrixDisplay = NULp;
66	static AP_userdef_matrix userdef_DNA_matrix(AP_MAX, AWAR_DIST_MATRIX_DNA_BASE);
67
68	static AP_matrix *get_user_matrix() {
69	AW_root *awr = AW_root::SINGLETON;
70	if (!awr->awar(AWAR_DIST_MATRIX_DNA_ENABLED)->read_int()) {
71	return NULp;
72	}
73	userdef_DNA_matrix.update_from_awars(awr);
74	userdef_DNA_matrix.normize();
75	return &userdef_DNA_matrix;
76	}
77
78	class BoundWindowCallback : virtual Noncopyable {
79	AW_window *aww;
80	WindowCallback cb;
81	public:
82	BoundWindowCallback(AW_window *aww_, const WindowCallback& cb_)
83	: aww(aww_),
84	cb(cb_)
85	{}
86	void operator()() { cb(aww); }
87	};
88
89	static SmartPtr<BoundWindowCallback> recalculate_matrix_cb;
90	static SmartPtr<BoundWindowCallback> recalculate_tree_cb;
91
92	static GB_ERROR last_matrix_calculation_error = NULp;
93
94	static void matrix_changed_cb() {
95	if (matrixDisplay) {
96	matrixDisplay->mark(MatrixDisplay::NEED_SETUP);
97	matrixDisplay->update_display();
98	}
99	}
100
101	struct RecalcNeeded {
102	bool matrix;
103	bool tree;
104
105	RecalcNeeded() : matrix(false), tree(false) { }
106	};
107
108	static RecalcNeeded need_recalc;
109
110	CONSTEXPR unsigned UPDATE_DELAY = 200;
111
112	static unsigned update_cb(AW_root *aw_root);
113	inline void add_update_cb() {
114	AW_root::SINGLETON->add_timed_callback(UPDATE_DELAY, makeTimedCallback(update_cb));
115	}
116
117	inline void matrix_needs_recalc_cb() {
118	need_recalc.matrix = true;
119	add_update_cb();
120	}
121	inline void tree_needs_recalc_cb() {
122	need_recalc.tree = true;
123	add_update_cb();
124	}
125
126	inline void compressed_matrix_needs_recalc_cb() {
127	if (GLOBAL_MATRIX.has_type(DI_MATRIX_COMPRESSED)) {
128	matrix_needs_recalc_cb();
129	}
130	}
131
132	static unsigned update_cb(AW_root *aw_root) {
133	if (need_recalc.matrix) {
134	GLOBAL_MATRIX.forget();
135	need_recalc.matrix = false; // because it's forgotten
136
137	int matrix_autocalc = aw_root->awar(AWAR_DIST_MATRIX_AUTO_RECALC)->read_int();
138	if (matrix_autocalc) {
139	bool recalc_now = true;
140	int tree_autocalc = aw_root->awar(AWAR_DIST_MATRIX_AUTO_CALC_TREE)->read_int();
141	if (!tree_autocalc) recalc_now = matrixDisplay ? matrixDisplay->willShow() : false;
142
143	if (recalc_now) {
144	di_assert(recalculate_matrix_cb.isSet());
145	(*recalculate_matrix_cb)();
146	di_assert(need_recalc.tree == true);
147	}
148	}
149	di_assert(need_recalc.matrix == false);
150	}
151
152	if (need_recalc.tree) {
153	int tree_autocalc = aw_root->awar(AWAR_DIST_MATRIX_AUTO_CALC_TREE)->read_int();
154	if (tree_autocalc) {
155	di_assert(recalculate_tree_cb.isSet());
156	(*recalculate_tree_cb)();
157	need_recalc.matrix = false; // otherwise endless loop, e.g. if output-tree is used for sorting
158	}
159	}
160
161	return 0; // do not call again
162	}
163
164	static void auto_calc_changed_cb(AW_root *aw_root) {
165	int matrix_autocalc = aw_root->awar(AWAR_DIST_MATRIX_AUTO_RECALC)->read_int();
166	int tree_autocalc = aw_root->awar(AWAR_DIST_MATRIX_AUTO_CALC_TREE)->read_int();
167
168	if (matrix_autocalc && !GLOBAL_MATRIX.exists()) matrix_needs_recalc_cb();
169	if (tree_autocalc && (matrix_autocalc \|\| GLOBAL_MATRIX.exists())) tree_needs_recalc_cb();
170	}
171
172	static AW_window create_dna_matrix_window(AW_root aw_root) {
173	AW_window_simple *aws = new AW_window_simple;
174	aws->init(aw_root, "SET_DNA_MATRIX", "SET MATRIX");
175	aws->auto_increment(50, 50);
176	aws->button_length(10);
177	aws->callback(AW_POPDOWN);
178	aws->create_button("CLOSE", "CLOSE");
179
180	aws->callback(makeHelpCallback("user_matrix.hlp"));
181	aws->create_button("HELP", "HELP");
182
183	aws->at_newline();
184
185	userdef_DNA_matrix.create_input_fields(aws);
186	aws->window_fit();
187	return aws;
188	}
189
190	static void selected_tree_changed_cb() {
191	if (AW_root::SINGLETON->awar(AWAR_DIST_CORR_TRANS)->read_int() == DI_TRANSFORMATION_FROM_TREE) {
192	matrix_needs_recalc_cb();
193	}
194	}
195
196	void DI_create_matrix_variables(AW_root *aw_root, AW_default def, AW_default db) {
197	GB_transaction ta(db);
198
199	userdef_DNA_matrix.set_descriptions(AP_A, "A");
200	userdef_DNA_matrix.set_descriptions(AP_C, "C");
201	userdef_DNA_matrix.set_descriptions(AP_G, "G");
202	userdef_DNA_matrix.set_descriptions(AP_T, "TU");
203	userdef_DNA_matrix.set_descriptions(AP_GAP, "GAP");
204
205	userdef_DNA_matrix.create_awars(aw_root);
206
207	RootCallback matrix_needs_recalc_callback = makeRootCallback(matrix_needs_recalc_cb);
208	aw_root->awar_int(AWAR_DIST_MATRIX_DNA_ENABLED, 0)->add_callback(matrix_needs_recalc_callback); // user matrix disabled by default
209	{
210	GBDATA *gbd = GB_search(AW_ROOT_DEFAULT, AWAR_DIST_MATRIX_DNA_BASE, GB_FIND);
211	GB_add_callback(gbd, GB_CB_CHANGED, makeDatabaseCallback(matrix_needs_recalc_cb));
212	}
213
214	aw_root->awar_string(AWAR_DIST_WHICH_SPECIES, "marked", def)->add_callback(matrix_needs_recalc_callback);
215	{
216	char *default_ali = GBT_get_default_alignment(db);
217	aw_root->awar_string(AWAR_DIST_ALIGNMENT, "", def)->add_callback(matrix_needs_recalc_callback)->write_string(default_ali);
218	free(default_ali);
219	}
220	aw_root->awar_string(AWAR_DIST_FILTER_ALIGNMENT, "none", def);
221	aw_root->awar_string(AWAR_DIST_FILTER_NAME, "none", def);
222	aw_root->awar_string(AWAR_DIST_FILTER_FILTER, "", def)->add_callback(matrix_needs_recalc_callback);
223	aw_root->awar_int (AWAR_DIST_FILTER_SIMPLIFY, 0, def)->add_callback(matrix_needs_recalc_callback);
224
225	aw_root->awar_string(AWAR_DIST_CANCEL_CHARS, ".", def)->add_callback(matrix_needs_recalc_callback);
226	aw_root->awar_int(AWAR_DIST_CORR_TRANS, (int)DI_TRANSFORMATION_SIMILARITY, def)->add_callback(matrix_needs_recalc_callback)->set_minmax(0, DI_TRANSFORMATION_COUNT-1);
227
228	aw_root->awar(AWAR_DIST_FILTER_ALIGNMENT)->map(AWAR_DIST_ALIGNMENT);
229
230	AW_create_fileselection_awars(aw_root, AWAR_DIST_SAVE_MATRIX_BASE, ".", "", "infile");
231	aw_root->awar_int(AWAR_DIST_SAVE_MATRIX_TYPE, 0, def);
232
233	enum treetype { CURR, SORT, COMPRESS, TREEAWARCOUNT };
234	AW_awar *tree_awar[TREEAWARCOUNT] = { NULp, NULp, NULp };
235
236	aw_root->awar_string(AWAR_DIST_TREE_STD_NAME, "tree_nj", def);
237	{
238	char *currentTree = aw_root->awar_string(AWAR_TREE, "", db)->read_string();
239	tree_awar[CURR] = aw_root->awar_string(AWAR_DIST_TREE_CURR_NAME, currentTree, def);
240	tree_awar[SORT] = aw_root->awar_string(AWAR_DIST_TREE_SORT_NAME, currentTree, def);
241	free(currentTree);
242	}
243	tree_awar[COMPRESS] = aw_root->awar_string(AWAR_DIST_TREE_COMP_NAME, NO_TREE_SELECTED, def);
244
245	aw_root->awar_int(AWAR_DIST_BOOTSTRAP_COUNT, 1000, def);
246	aw_root->awar_int(AWAR_DIST_MATRIX_AUTO_RECALC, 0, def)->add_callback(auto_calc_changed_cb);
247	aw_root->awar_int(AWAR_DIST_MATRIX_AUTO_CALC_TREE, 0, def)->add_callback(auto_calc_changed_cb);
248
249	aw_root->awar_float(AWAR_DIST_MIN_DIST, 0.0)->set_minmax(0.0, 4.0);
250	aw_root->awar_float(AWAR_DIST_MAX_DIST, 0.0)->set_minmax(0.0, 4.0);
251
252	aw_root->awar_string(AWAR_SPECIES_NAME, "", db);
253
254	DI_create_cluster_awars(aw_root, def, db);
255
256	#if defined(DEBUG)
257	AWT_create_db_browser_awars(aw_root, def);
258	#endif // DEBUG
259
260	{
261	GB_push_transaction(db);
262
263	GBDATA *gb_species_data = GBT_get_species_data(db);
264	GB_add_callback(gb_species_data, GB_CB_CHANGED, makeDatabaseCallback(matrix_needs_recalc_cb));
265
266	GB_pop_transaction(db);
267	}
268
269	AWT_registerTreeAwarCallback(tree_awar[CURR], makeTreeAwarCallback(selected_tree_changed_cb), true);
270	AWT_registerTreeAwarCallback(tree_awar[SORT], makeTreeAwarCallback(matrix_needs_recalc_cb), true);
271	AWT_registerTreeAwarCallback(tree_awar[COMPRESS], makeTreeAwarCallback(compressed_matrix_needs_recalc_cb), true);
272
273	auto_calc_changed_cb(aw_root);
274	}
275
276	DI_ENTRY::DI_ENTRY(GBDATA gbd, DI_MATRIX phmatrix_)
277	: phmatrix(phmatrix_),
278	full_name(NULp),
279	sequence(NULp),
280	name(NULp),
281	group_nr(0)
282	{
283	GBDATA *gb_ali = GB_entry(gbd, phmatrix->get_aliname());
284	if (gb_ali) {
285	GBDATA *gb_data = GB_entry(gb_ali, "data");
286	if (gb_data) {
287	if (phmatrix->is_AA) {
288	sequence = new AP_sequence_simple_protein(phmatrix->get_aliview());
289	}
290	else {
291	sequence = new AP_sequence_parsimony(phmatrix->get_aliview());
292	}
293	sequence->bind_to_species(gbd);
294	sequence->lazy_load_sequence(); // load sequence
295
296	name = GBT_read_string(gbd, "name");
297	full_name = GBT_read_string(gbd, "full_name");
298	}
299	}
300	}
301
302	DI_ENTRY::DI_ENTRY(const char name_, DI_MATRIX phmatrix_)
303	: phmatrix(phmatrix_),
304	full_name(NULp),
305	sequence(NULp),
306	name(ARB_strdup(name_)),
307	group_nr(0)
308	{}
309
310	DI_ENTRY::~DI_ENTRY() {
311	delete sequence;
312	free(name);
313	free(full_name);
314	}
315
316	DI_MATRIX::DI_MATRIX(const AliView& aliview_)
317	: gb_species_data(NULp),
318	seq_len(0),
319	allocated_entries(0),
320	aliview(new AliView(aliview_)),
321	is_AA(false),
322	entries(NULp),
323	nentries(0),
324	matrix(NULp),
325	matrix_type(DI_MATRIX_FULL)
326	{
327	memset(cancel_columns, 0, sizeof(cancel_columns));
328	}
329
330	char *DI_MATRIX::unload() {
331	for (size_t i=0; i<nentries; i++) {
332	delete entries[i];
333	}
334	freenull(entries);
335	nentries = 0;
336	return NULp;
337	}
338
339	DI_MATRIX::~DI_MATRIX() {
340	unload();
341	delete matrix;
342	delete aliview;
343	}
344
345	struct TreeOrderedSpecies {
346	GBDATA *gbd;
347	int order_index;
348
349	TreeOrderedSpecies(const MatrixOrder& order, GBDATA *gb_spec)
350	: gbd(gb_spec),
351	order_index(order.get_index(GBT_get_name_or_description(gbd)))
352	{}
353	};
354
355	MatrixOrder::MatrixOrder(GBDATA *gb_main, GB_CSTR sort_tree_name)
356	: name2pos(NULp),
357	leafs(0)
358	{
359	if (sort_tree_name) {
360	int size;
361	TreeNode *sort_tree = GBT_read_tree_and_size(gb_main, sort_tree_name, new SimpleRoot, &size);
362
363	if (sort_tree) {
364	leafs = size+1;
365	name2pos = GBS_create_hash(leafs, GB_IGNORE_CASE);
366
367	IF_ASSERTION_USED(int leafsLoaded = leafs);
368	leafs = 0;
369	insert_in_hash(sort_tree);
370
371	arb_assert(leafsLoaded == leafs);
372	destroy(sort_tree);
373	}
374	else {
375	GB_clear_error();
376	}
377	}
378	}
379	static int TreeOrderedSpecies_cmp(const void p1, const void p2, void *) {
380	TreeOrderedSpecies s1 = (TreeOrderedSpecies)p1;
381	TreeOrderedSpecies s2 = (TreeOrderedSpecies)p2;
382
383	return s2->order_index - s1->order_index;
384	}
385
386	void MatrixOrder::applyTo(TreeOrderedSpecies **species_array, size_t array_size) const {
387	GB_sort((void**)species_array, 0, array_size, TreeOrderedSpecies_cmp, NULp);
388	}
389
390	GB_ERROR DI_MATRIX::load(LoadWhat what, const MatrixOrder& order, bool show_warnings, GBDATA **species_list) {
391	GBDATA *gb_main = get_gb_main();
392	const char *use = get_aliname();
393
394	GB_transaction ta(gb_main);
395
396	seq_len = GBT_get_alignment_len(gb_main, use);
397	is_AA = GBT_is_alignment_protein(gb_main, use);
398	gb_species_data = GBT_get_species_data(gb_main);
399
400	allocated_entries = 1000;
401	ARB_calloc(entries, allocated_entries);
402
403	nentries = 0;
404
405	size_t no_of_species = -1U;
406	switch (what) {
407	case DI_LOAD_ALL:
408	no_of_species = GBT_get_species_count(gb_main);
409	break;
410	case DI_LOAD_MARKED:
411	no_of_species = GBT_count_marked_species(gb_main);
412	break;
413	case DI_LOAD_LIST:
414	di_assert(species_list);
415	for (no_of_species = 0; species_list[no_of_species]; ++no_of_species) ;
416	break;
417	}
418
419	di_assert(no_of_species != -1U);
420	if (no_of_species<2) {
421	return GBS_global_string("Not enough input species (%zu)", no_of_species);
422	}
423
424	TreeOrderedSpecies *species_to_load[no_of_species];
425
426	{
427	size_t i = 0;
428	switch (what) {
429	case DI_LOAD_ALL: {
430	for (GBDATA *gb_species = GBT_first_species_rel_species_data(gb_species_data); gb_species; gb_species = GBT_next_species(gb_species), ++i) {
431	species_to_load[i] = new TreeOrderedSpecies(order, gb_species);
432	}
433	break;
434	}
435	case DI_LOAD_MARKED: {
436	for (GBDATA *gb_species = GBT_first_marked_species_rel_species_data(gb_species_data); gb_species; gb_species = GBT_next_marked_species(gb_species), ++i) {
437	species_to_load[i] = new TreeOrderedSpecies(order, gb_species);
438	}
439	break;
440	}
441	case DI_LOAD_LIST: {
442	for (i = 0; species_list[i]; ++i) {
443	species_to_load[i] = new TreeOrderedSpecies(order, species_list[i]);
444	}
445	break;
446	}
447	}
448	arb_assert(i == no_of_species);
449	}
450
451	if (order.defined()) {
452	order.applyTo(species_to_load, no_of_species);
453	if (show_warnings) {
454	int species_not_in_sort_tree = 0;
455	for (size_t i = 0; i<no_of_species; ++i) {
456	if (!species_to_load[i]->order_index) {
457	species_not_in_sort_tree++;
458	}
459	}
460	if (species_not_in_sort_tree) {
461	GBT_message(gb_main, GBS_global_string("Warning: %i of the affected species are not in sort-tree", species_not_in_sort_tree));
462	}
463	}
464	}
465	else {
466	if (show_warnings) {
467	static bool shown = false;
468	if (!shown) { // showing once is enough
469	GBT_message(gb_main, "Warning: No valid tree given to sort matrix (using default database order)");
470	shown = true;
471	}
472	}
473	}
474
475	if (no_of_species>allocated_entries) {
476	allocated_entries = no_of_species;
477	ARB_realloc(entries, allocated_entries);
478	}
479
480	GB_ERROR error = NULp;
481	arb_progress progress("Preparing sequence data", no_of_species);
482	for (size_t i = 0; i<no_of_species && !error; ++i) {
483	DI_ENTRY *phentry = new DI_ENTRY(species_to_load[i]->gbd, this);
484	if (phentry->sequence) { // a species found
485	arb_assert(nentries<allocated_entries);
486	entries[nentries++] = phentry;
487	}
488	else {
489	delete phentry;
490	}
491	delete species_to_load[i];
492	species_to_load[i] = NULp;
493
494	progress.inc_and_check_user_abort(error);
495	}
496
497	return error;
498	}
499
500	char DI_MATRIX::calculate_overall_freqs(double rel_frequencies[AP_MAX], char cancel) {
501	di_assert(is_AA == false);
502
503	long hits2[AP_MAX];
504	long sum = 0;
505	int i;
506	int pos;
507	int b;
508	long s_len = aliview->get_length();
509
510	memset((char *) &hits2[0], 0, sizeof(hits2));
511	for (size_t row = 0; row < nentries; row++) {
512	const char *seq1 = entries[row]->get_nucl_seq()->get_sequence();
513	for (pos = 0; pos < s_len; pos++) {
514	b = *(seq1++);
515	if (cancel[b]) continue;
516	hits2[b]++;
517	}
518	}
519	for (i = 0; i < AP_MAX; i++) { // LOOP_VECTORIZED
520	sum += hits2[i];
521	}
522	for (i = 0; i < AP_MAX; i++) {
523	rel_frequencies[i] = hits2[i] / (double) sum;
524	}
525	return NULp;
526	}
527
528	double DI_MATRIX::corr(double dist, double b, double & sigma) {
529	const double eps = 0.01;
530	double ar = 1.0 - dist/b;
531	sigma = 1000.0;
532	if (ar< eps) return 3.0;
533	sigma = b/ar;
534	return - b * log(1-dist/b);
535	}
536
537	GB_ERROR DI_MATRIX::calculate(const char cancel, DI_TRANSFORMATION transformation, bool aborted_flag, AP_matrix *userdef_matrix) {
538	di_assert(is_AA == false);
539
540	if (userdef_matrix) {
541	switch (transformation) {
542	case DI_TRANSFORMATION_NONE:
543	case DI_TRANSFORMATION_SIMILARITY:
544	case DI_TRANSFORMATION_JUKES_CANTOR:
545	break;
546	default:
547	aw_message("Sorry: this kind of distance correction does not support a user defined matrix - it will be ignored");
548	userdef_matrix = NULp;
549	break;
550	}
551	}
552
553	matrix = new AP_smatrix(nentries);
554
555	long s_len = aliview->get_length();
556	long hits[AP_MAX][AP_MAX];
557	size_t i;
558
559	if (nentries<=1) {
560	return "Not enough species selected to calculate matrix";
561	}
562	memset(&cancel_columns[0], 0, 256);
563
564	for (i=0; cancel[i]; i++) {
565	cancel_columns[safeCharIndex(AP_sequence_parsimony::table[safeCharIndex(cancel[i])])] = 1;
566	UNCOVERED(); // @@@ cover
567	}
568
569	long columns;
570	double b;
571	long frequencies[AP_MAX];
572	double rel_frequencies[AP_MAX];
573	double S_square = 0;
574
575	switch (transformation) {
576	case DI_TRANSFORMATION_FELSENSTEIN:
577	this->calculate_overall_freqs(rel_frequencies, cancel_columns);
578	S_square = 0.0;
579	for (i=0; i<AP_MAX; i++) { // LOOP_VECTORIZED[!<8.1]
580	S_square += rel_frequencies[i]*rel_frequencies[i];
581	}
582	break;
583	default: break;
584	};
585
586	arb_progress progress("Calculating distance matrix", matrix_halfsize(nentries, true));
587	GB_ERROR error = NULp;
588	for (size_t row = 0; row<nentries && !error; row++) {
589	for (size_t col=0; col<=row && !error; col++) {
590	columns = 0;
591
592	const unsigned char *seq1 = entries[row]->get_nucl_seq()->get_usequence();
593	const unsigned char *seq2 = entries[col]->get_nucl_seq()->get_usequence();
594
595	b = 0.0;
596	switch (transformation) {
597	case DI_TRANSFORMATION_FROM_TREE:
598	di_assert(0);
599	break;
600	case DI_TRANSFORMATION_JUKES_CANTOR:
601	b = 0.75;
602	// fall-through
603	case DI_TRANSFORMATION_NONE:
604	case DI_TRANSFORMATION_SIMILARITY: {
605	double dist = 0.0;
606	if (userdef_matrix) {
607	memset((char )hits, 0, sizeof(long) AP_MAX * AP_MAX);
608	int pos;
609	for (pos = s_len; pos >= 0; pos--) {
610	hits[(seq1++)][(seq2++)]++;
611	}
612	int x, y;
613	double diffsum = 0.0;
614	double all_sum = 0.001;
615	for (x = AP_A; x < AP_MAX; x*=2) {
616	for (y = AP_A; y < AP_MAX; y*=2) {
617	if (x==y) {
618	all_sum += hits[x][y];
619	}
620	else {
621	UNCOVERED(); // @@@ cover
622	diffsum += hits[x][y] * userdef_matrix->get(x, y);
623	all_sum += hits[x][y] * userdef_matrix->get(x, y);
624	}
625	}
626	}
627	dist = diffsum / all_sum;
628	}
629	else {
630	for (int pos = s_len; pos >= 0; pos--) {
631	int b1 = *(seq1++);
632	int b2 = *(seq2++);
633	if (cancel_columns[b1]) continue;
634	if (cancel_columns[b2]) continue;
635	columns++;
636	if (b1&b2) continue;
637	dist+=1.0;
638	}
639	if (columns == 0) columns = 1;
640	dist /= columns;
641	}
642	if (transformation==DI_TRANSFORMATION_SIMILARITY) {
643	dist = (1.0-dist);
644	}
645	else if (b) {
646	double sigma;
647	dist = this->corr(dist, b, sigma);
648	}
649	matrix->set(row, col, dist);
650	break;
651	}
652	case DI_TRANSFORMATION_KIMURA:
653	case DI_TRANSFORMATION_OLSEN:
654	case DI_TRANSFORMATION_FELSENSTEIN: {
655	int pos;
656	double dist = 0.0;
657	long N, P, Q, M;
658	double p, q;
659
660	memset((char )hits, 0, sizeof(long) AP_MAX * AP_MAX);
661	for (pos = s_len; pos >= 0; pos--) {
662	hits[(seq1++)][(seq2++)]++;
663	}
664	switch (transformation) {
665	case DI_TRANSFORMATION_KIMURA:
666	P = hits[AP_A][AP_G] +
667	hits[AP_G][AP_A] +
668	hits[AP_C][AP_T] +
669	hits[AP_T][AP_C];
670	Q = hits[AP_A][AP_C] +
671	hits[AP_A][AP_T] +
672	hits[AP_C][AP_A] +
673	hits[AP_T][AP_A] +
674	hits[AP_G][AP_C] +
675	hits[AP_G][AP_T] +
676	hits[AP_C][AP_G] +
677	hits[AP_T][AP_G];
678	M = hits[AP_A][AP_A] +
679	hits[AP_C][AP_C] +
680	hits[AP_G][AP_G] +
681	hits[AP_T][AP_T];
682	N = P+Q+M;
683	if (N==0) N=1;
684	p = (double)P/(double)N;
685	q = (double)Q/(double)N;
686	dist = - .5 * log(
687	(1.0-2.0p-q)sqrt(1.0-2.0*q)
688	);
689	break;
690
691	case DI_TRANSFORMATION_OLSEN:
692	case DI_TRANSFORMATION_FELSENSTEIN:
693
694	memset((char *)frequencies, 0,
695	sizeof(long) * AP_MAX);
696
697	N = 0;
698	M = 0;
699
700	for (i=0; i<AP_MAX; i++) {
701	if (cancel_columns[i]) continue;
702	unsigned int j;
703	for (j=0; j<i; j++) {
704	if (cancel_columns[j]) continue;
705	frequencies[i] +=
706	hits[i][j]+
707	hits[j][i];
708	}
709	frequencies[i] += hits[i][i];
710	N += frequencies[i];
711	M += hits[i][i];
712	}
713	if (N==0) N=1;
714	if (transformation == DI_TRANSFORMATION_OLSEN) { // Calc sum square freq individually for each line
715	S_square = 0.0;
716	for (i=0; i<AP_MAX; i++) S_square += frequencies[i]*frequencies[i];
717	b = 1.0 - S_square/((double)N*(double)N);
718	}
719	else {
720	b = 1.0 - S_square;
721	}
722
723	dist = ((double)(N-M)) / (double) N;
724	double sigma;
725	dist = this->corr(dist, b, sigma);
726	break;
727
728	default: return "Sorry: Transformation not implemented";
729	}
730	matrix->set(row, col, dist);
731	break;
732	}
733	default:;
734	} // switch
735	progress.inc_and_check_user_abort(error);
736	}
737	}
738	if (aborted_flag && progress.aborted()) *aborted_flag = true;
739	return error;
740	}
741
742	GB_ERROR DI_MATRIX::calculate_pro(DI_TRANSFORMATION transformation, bool *aborted_flag) {
743	di_assert(is_AA == true);
744
745	di_cattype catType;
746	switch (transformation) {
747	case DI_TRANSFORMATION_NONE: catType = NONE; break;
748	case DI_TRANSFORMATION_SIMILARITY: catType = SIMILARITY; break;
749	case DI_TRANSFORMATION_KIMURA: catType = KIMURA; break;
750	case DI_TRANSFORMATION_PAM: catType = PAM; break;
751	case DI_TRANSFORMATION_CATEGORIES_HALL: catType = HALL; break;
752	case DI_TRANSFORMATION_CATEGORIES_BARKER: catType = GEORGE; break;
753	case DI_TRANSFORMATION_CATEGORIES_CHEMICAL: catType = CHEMICAL; break;
754	default:
755	return "This correction is not available for protein data";
756	}
757	matrix = new AP_smatrix(nentries);
758
759	di_protdist prodist(UNIVERSAL, catType, nentries, entries, aliview->get_length(), matrix);
760	return prodist.makedists(aborted_flag);
761	}
762
763	typedef std::map<const char, TreeNode, charpLess> NamedNodes;
764
765	GB_ERROR link_to_tree(NamedNodes& named, TreeNode *node) {
766	GB_ERROR error = NULp;
767	if (node->is_leaf()) {
768	NamedNodes::iterator found = named.find(node->name);
769	if (found != named.end()) {
770	if (found->second) {
771	error = GBS_global_string("Invalid tree (two nodes named '%s')", node->name);
772	}
773	else {
774	found->second = node;
775	}
776	}
777	// otherwise, we do not care about the node (e.g. because it is not marked)
778	}
779	else {
780	error = link_to_tree(named, node->get_leftson());
781	if (!error) error = link_to_tree(named, node->get_rightson());
782	}
783	return error;
784	}
785
786	#if defined(ASSERTION_USED)
787	static TreeNode findNode(TreeNode node, const char *name) {
788	if (node->is_leaf()) {
789	return strcmp(node->name, name) == 0 ? node : NULp;
790	}
791
792	TreeNode *found = findNode(node->get_leftson(), name);
793	if (!found) found = findNode(node->get_rightson(), name);
794	return found;
795	}
796	#endif
797
798	static GB_ERROR init(NamedNodes& node, TreeNode tree, const DI_ENTRYconst*const entries, size_t nentries) {
799	GB_ERROR error = NULp;
800	for (size_t n = 0; n<nentries; ++n) {
801	node[entries[n]->name] = NULp;
802	}
803	error = link_to_tree(node, tree);
804	if (!error) { // check for missing species (needed but not in tree)
805	size_t missing = 0;
806	const char *exampleName = NULp;
807
808	for (size_t n = 0; n<nentries; ++n) {
809	NamedNodes::iterator found = node.find(entries[n]->name);
810	if (found == node.end()) {
811	++missing;
812	exampleName = entries[n]->name;
813	}
814	else {
815	di_assert(node[entries[n]->name] == findNode(tree, entries[n]->name));
816	if (!node[entries[n]->name]) {
817	++missing;
818	exampleName = entries[n]->name;
819	}
820	}
821	}
822
823	if (missing) {
824	error = GBS_global_string("Tree is missing %zu required species (e.g. '%s')", missing, exampleName);
825	}
826	}
827	return error;
828	}
829
830	GB_ERROR DI_MATRIX::extract_from_tree(const char treename, bool aborted_flag) {
831	GB_ERROR error = NULp;
832	if (nentries<=1) error = "Not enough species selected to calculate matrix";
833	else {
834	TreeNode *tree;
835	{
836	GB_transaction ta(get_gb_main());
837	tree = GBT_read_tree(get_gb_main(), treename, new SimpleRoot);
838	}
839	if (!tree) error = GB_await_error();
840	else {
841	arb_progress progress("Extracting distances from tree", matrix_halfsize(nentries, true));
842	NamedNodes node;
843
844	error = init(node, tree, entries, nentries);
845	matrix = new AP_smatrix(nentries);
846
847	for (size_t row = 0; row<nentries && !error; row++) {
848	TreeNode *rnode = node[entries[row]->name];
849	for (size_t col=0; col<=row && !error; col++) {
850	double dist;
851	if (col != row) {
852	TreeNode *cnode = node[entries[col]->name];
853	dist = rnode->intree_distance_to(cnode);
854	}
855	else {
856	dist = 0.0;
857	}
858	matrix->set(row, col, dist);
859	progress.inc_and_check_user_abort(error);
860	}
861	}
862	UNCOVERED();
863	destroy(tree);
864	if (aborted_flag && progress.aborted()) *aborted_flag = true;
865	if (error) progress.done();
866	}
867	}
868	return error;
869	}
870
871	__ATTR__USERESULT static GB_ERROR di_calculate_matrix(AW_root aw_root, const WeightedFilter weighted_filter, bool bootstrap_flag, bool show_warnings, bool *aborted_flag) {
872	// sets 'aborted_flag' to true, if it is non-NULp and the calculation has been aborted
873	GB_ERROR error = NULp;
874
875	if (GLOBAL_MATRIX.exists()) {
876	di_assert(!need_recalc.matrix);
877	}
878	else {
879	GBDATA *gb_main = weighted_filter->get_gb_main();
880	GB_transaction ta(gb_main);
881
882	char *use = aw_root->awar(AWAR_DIST_ALIGNMENT)->read_string();
883	long ali_len = GBT_get_alignment_len(gb_main, use);
884
885	if (ali_len<=0) {
886	error = "Please select a valid alignment";
887	GB_clear_error();
888	}
889	else {
890	arb_progress progress("Calculating matrix");
891	char *cancel = aw_root->awar(AWAR_DIST_CANCEL_CHARS)->read_string();
892	AliView *aliview = weighted_filter->create_aliview(use, error);
893
894	if (!error) {
895	if (bootstrap_flag) aliview->get_filter()->enable_bootstrap();
896
897	char *load_what = aw_root->awar(AWAR_DIST_WHICH_SPECIES)->read_string();
898	char *sort_tree_name = aw_root->awar(AWAR_DIST_TREE_SORT_NAME)->read_string();
899
900	LoadWhat all_flag = (strcmp(load_what, "all") == 0) ? DI_LOAD_ALL : DI_LOAD_MARKED;
901	{
902	DI_MATRIX phm = new DI_MATRIX(aliview);
903	phm->matrix_type = DI_MATRIX_FULL;
904
905	static SmartCharPtr last_sort_tree_name;
906	static SmartPtr<MatrixOrder> last_order;
907
908	if (last_sort_tree_name.isNull() \|\| !sort_tree_name \|\| strcmp(&*last_sort_tree_name, sort_tree_name) != 0) {
909	last_sort_tree_name = nulldup(sort_tree_name);
910	last_order = new MatrixOrder(gb_main, sort_tree_name);
911	}
912	di_assert(last_order.isSet());
913	error = phm->load(all_flag, *last_order, show_warnings, NULp);
914
915	free(sort_tree_name);
916	error = ta.close(error);
917
918	bool aborted = false;
919	if (!error) {
920	if (progress.aborted()) {
921	phm->unload();
922	error = "Aborted by user";
923	aborted = true;
924	}
925	else {
926	DI_TRANSFORMATION trans = (DI_TRANSFORMATION)aw_root->awar(AWAR_DIST_CORR_TRANS)->read_int();
927
928	if (trans == DI_TRANSFORMATION_FROM_TREE) {
929	const char *treename = aw_root->awar(AWAR_DIST_TREE_CURR_NAME)->read_char_pntr();
930	error = phm->extract_from_tree(treename, &aborted);
931	}
932	else {
933	if (phm->is_AA) error = phm->calculate_pro(trans, &aborted);
934	else error = phm->calculate(cancel, trans, &aborted, get_user_matrix());
935	}
936	}
937	}
938
939	if (aborted) {
940	di_assert(error);
941	if (aborted_flag) *aborted_flag = true;
942	}
943	if (error) {
944	delete phm;
945	GLOBAL_MATRIX.forget();
946	}
947	else {
948	GLOBAL_MATRIX.replaceBy(phm);
949	tree_needs_recalc_cb();
950	need_recalc.matrix = false;
951	}
952	}
953	free(load_what);
954	}
955
956	free(cancel);
957	delete aliview;
958	}
959	free(use);
960
961	di_assert(contradicted(error, GLOBAL_MATRIX.exists()));
962	}
963	return error;
964	}
965
966	static void di_mark_by_distance(AW_window aww, WeightedFilter weighted_filter) {
967	AW_root *aw_root = aww->get_root();
968	float lowerBound = aw_root->awar(AWAR_DIST_MIN_DIST)->read_float();
969	float upperBound = aw_root->awar(AWAR_DIST_MAX_DIST)->read_float();
970
971	GB_ERROR error = NULp;
972	if (lowerBound >= upperBound) {
973	error = GBS_global_string("Lower bound (%f) has to be smaller than upper bound (%f)", lowerBound, upperBound);
974	}
975	else if (lowerBound<0.0 \|\| lowerBound > 1.0) {
976	error = GBS_global_string("Lower bound (%f) is not in allowed range [0.0 .. 1.0]", lowerBound);
977	}
978	else if (upperBound<0.0 \|\| upperBound > 1.0) {
979	error = GBS_global_string("Upper bound (%f) is not in allowed range [0.0 .. 1.0]", upperBound);
980	}
981	else {
982	GBDATA *gb_main = weighted_filter->get_gb_main();
983	GB_transaction ta(gb_main);
984
985	char *selected = aw_root->awar(AWAR_SPECIES_NAME)->read_string();
986	if (!selected[0]) {
987	error = "Please select a species";
988	}
989	else {
990	GBDATA *gb_selected = GBT_find_species(gb_main, selected);
991	if (!gb_selected) {
992	error = GBS_global_string("Couldn't find species '%s'", selected);
993	}
994	else {
995	char *use = aw_root->awar(AWAR_DIST_ALIGNMENT)->read_string();
996	char *cancel = aw_root->awar(AWAR_DIST_CANCEL_CHARS)->read_string();
997	DI_TRANSFORMATION trans = (DI_TRANSFORMATION)aw_root->awar(AWAR_DIST_CORR_TRANS)->read_int();
998	AliView *aliview = weighted_filter->create_aliview(use, error);
999
1000	if (!error) {
1001	DI_MATRIX *prev_global = GLOBAL_MATRIX.swap(NULp);
1002
1003	size_t speciesCount = GBT_get_species_count(gb_main);
1004	bool markedSelected = false;
1005
1006	arb_progress progress("Mark species by distance", speciesCount);
1007	MatrixOrder order(gb_main, NULp);
1008
1009	for (GBDATA *gb_species = GBT_first_species(gb_main);
1010	gb_species && !error;
1011	gb_species = GBT_next_species(gb_species))
1012	{
1013	DI_MATRIX phm = new DI_MATRIX(aliview);
1014	phm->matrix_type = DI_MATRIX_FULL;
1015	GBDATA *species_pair[] = { gb_selected, gb_species, NULp };
1016
1017	error = phm->load(DI_LOAD_LIST, order, false, species_pair);
1018
1019	if (phm->nentries == 2) { // if species has no alignment -> nentries<2
1020	if (!error) {
1021	if (phm->is_AA) error = phm->calculate_pro(trans, NULp);
1022	else error = phm->calculate(cancel, trans, NULp, get_user_matrix());
1023	}
1024
1025	if (!error) {
1026	double dist_value = phm->matrix->get(0, 1); // distance or conformance
1027	bool mark = (lowerBound <= dist_value && dist_value <= upperBound);
1028	GB_write_flag(gb_species, mark);
1029
1030	if (!markedSelected) {
1031	dist_value = phm->matrix->get(0, 0); // distance or conformance to self
1032	mark = (lowerBound <= dist_value && dist_value <= upperBound);
1033	GB_write_flag(gb_selected, mark);
1034
1035	markedSelected = true;
1036	}
1037	}
1038	}
1039
1040	delete phm;
1041	if (!error) progress.inc_and_check_user_abort(error);
1042	}
1043
1044	di_assert(!GLOBAL_MATRIX.exists());
1045	ASSERT_RESULT(DI_MATRIX*, NULp, GLOBAL_MATRIX.swap(prev_global));
1046
1047	if (error) progress.done();
1048	}
1049
1050	delete aliview;
1051	free(cancel);
1052	free(use);
1053	}
1054	}
1055
1056	free(selected);
1057	error = ta.close(error);
1058	}
1059
1060	if (error) {
1061	aw_message(error);
1062	}
1063	}
1064
1065	static GB_ERROR di_recalc_matrix() {
1066	// recalculate matrix
1067	last_matrix_calculation_error = NULp;
1068	if (need_recalc.matrix && GLOBAL_MATRIX.exists()) {
1069	GLOBAL_MATRIX.forget();
1070	}
1071	di_assert(recalculate_matrix_cb.isSet());
1072	(*recalculate_matrix_cb)();
1073	return last_matrix_calculation_error;
1074	}
1075
1076	static void di_view_matrix_cb(AW_window aww, save_matrix_params sparam) {
1077	GB_ERROR error = di_recalc_matrix();
1078	if (!error) {
1079	if (!matrixDisplay) matrixDisplay = new MatrixDisplay(sparam->weighted_filter->get_gb_main());
1080
1081	static AW_window *viewer = NULp;
1082	if (!viewer) viewer = DI_create_view_matrix_window(aww->get_root(), matrixDisplay, sparam);
1083
1084	matrixDisplay->mark(MatrixDisplay::NEED_SETUP);
1085	matrixDisplay->update_display();
1086
1087	GLOBAL_MATRIX.set_changed_cb(matrix_changed_cb);
1088
1089	viewer->activate();
1090	}
1091	}
1092
1093	static void di_save_matrix_cb(AW_window *aww) {
1094	// save the matrix
1095	GB_ERROR error = di_recalc_matrix();
1096	if (!error) {
1097	char *filename = aww->get_root()->awar(AWAR_DIST_SAVE_MATRIX_FILENAME)->read_string();
1098	enum DI_SAVE_TYPE type = (enum DI_SAVE_TYPE)aww->get_root()->awar(AWAR_DIST_SAVE_MATRIX_TYPE)->read_int();
1099
1100	GLOBAL_MATRIX.get()->save(filename, type);
1101	free(filename);
1102	}
1103	AW_refresh_fileselection(aww->get_root(), AWAR_DIST_SAVE_MATRIX_BASE);
1104	aww->hide_or_notify(error);
1105	}
1106
1107	AW_window DI_create_save_matrix_window(AW_root aw_root, save_matrix_params *save_params) {
1108	static AW_window_simple *aws = NULp;
1109	if (!aws) {
1110	aws = new AW_window_simple;
1111	aws->init(aw_root, "SAVE_MATRIX", "Save Matrix");
1112	aws->load_xfig("sel_box_user.fig");
1113
1114	aws->at("close");
1115	aws->callback(AW_POPDOWN);
1116	aws->create_button("CLOSE", "CANCEL", "C");
1117
1118
1119	aws->at("help"); aws->callback(makeHelpCallback("save_matrix.hlp"));
1120	aws->create_button("HELP", "HELP", "H");
1121
1122	aws->at("user");
1123	aws->create_option_menu(AWAR_DIST_SAVE_MATRIX_TYPE, true);
1124	aws->insert_default_option("Phylip Format (Lower Triangular Matrix)", "P", DI_SAVE_PHYLIP_COMP);
1125	aws->insert_option("Readable (using NDS)", "R", DI_SAVE_READABLE);
1126	aws->insert_option("Tabbed (using NDS)", "R", DI_SAVE_TABBED);
1127	aws->update_option_menu();
1128
1129	AW_create_standard_fileselection(aws, save_params->awar_base);
1130
1131	aws->at("save2");
1132	aws->callback(makeWindowCallback(di_save_matrix_cb));
1133	aws->create_button("SAVE", "SAVE", "S");
1134
1135	aws->at("cancel2");
1136	aws->callback(AW_POPDOWN);
1137	aws->create_button("CLOSE", "CANCEL", "C");
1138	}
1139	return aws;
1140	}
1141
1142	static AW_window awt_create_select_cancel_window(AW_root aw_root) {
1143	AW_window_simple *aws = new AW_window_simple;
1144	aws->init(aw_root, "SELECT_CHARS_TO_CANCEL_COLUMN", "CANCEL SELECT");
1145	aws->load_xfig("di_cancel.fig");
1146
1147	aws->at("close");
1148	aws->callback(AW_POPDOWN);
1149	aws->create_button("CLOSE", "CLOSE", "C");
1150
1151	aws->at("cancel");
1152	aws->create_input_field(AWAR_DIST_CANCEL_CHARS, 12);
1153
1154	return aws;
1155	}
1156
1157	static const char *enum_trans_to_string[] = {
1158	"none",
1159	"similarity",
1160	"jukes_cantor",
1161	"felsenstein",
1162
1163	"pam",
1164	"hall",
1165	"barker",
1166	"chemical",
1167
1168	"kimura",
1169	"olsen",
1170	"felsenstein voigt",
1171	"olsen voigt",
1172	"max ml",
1173
1174	NULp, // treedist
1175	};
1176
1177	STATIC_ASSERT(ARRAY_ELEMS(enum_trans_to_string) == DI_TRANSFORMATION_COUNT);
1178
1179	static void di_calculate_tree_cb(AW_window aww, WeightedFilter weighted_filter, bool bootstrap_flag) {
1180	recalculate_tree_cb = new BoundWindowCallback(aww, makeWindowCallback(di_calculate_tree_cb, weighted_filter, bootstrap_flag));
1181
1182	AW_root *aw_root = aww->get_root();
1183	GB_ERROR error = NULp;
1184	StrArray *all_names = NULp;
1185
1186	long loop_count = 0;
1187	long bootstrap_count = aw_root->awar(AWAR_DIST_BOOTSTRAP_COUNT)->read_int();
1188
1189	{
1190	char *tree_name = aw_root->awar(AWAR_DIST_TREE_STD_NAME)->read_string();
1191	error = GBT_check_tree_name(tree_name);
1192	free(tree_name);
1193	}
1194
1195	SmartPtr<arb_progress> progress;
1196	SmartPtr<ConsensusTree> ctree;
1197
1198	if (!error) {
1199	if (bootstrap_flag) {
1200	if (bootstrap_count) {
1201	progress = new arb_progress("Calculating bootstrap trees", bootstrap_count+1);
1202	}
1203	else {
1204	progress = new arb_progress("Calculating bootstrap trees (KILL to stop)", long(INT_MAX));
1205	}
1206	progress->auto_subtitles("tree");
1207	}
1208	else {
1209	progress = new arb_progress("Calculating tree");
1210	}
1211
1212	if (bootstrap_flag) {
1213	GLOBAL_MATRIX.forget();
1214	GLOBAL_MATRIX.set_changed_cb(NULp); // otherwise matrix window will repeatedly pop up/down
1215
1216	error = di_calculate_matrix(aw_root, weighted_filter, bootstrap_flag, true, NULp);
1217	if (!error) {
1218	DI_MATRIX *matr = GLOBAL_MATRIX.get();
1219	if (!matr) {
1220	error = "unexpected error in di_calculate_matrix_cb (data missing)";
1221	}
1222	else {
1223	all_names = new StrArray;
1224	all_names->reserve(matr->nentries+2);
1225
1226	for (size_t i=0; i<matr->nentries; i++) {
1227	all_names->put(ARB_strdup(matr->entries[i]->name));
1228	}
1229	ctree = new ConsensusTree(*all_names);
1230	}
1231	}
1232	}
1233	}
1234
1235	TreeNode *tree = NULp;
1236	do {
1237	if (error) break;
1238
1239	bool aborted = false;
1240
1241	if (bootstrap_flag) {
1242	if (loop_count>0) { // in first loop we already have a valid matrix -> no need to recalculate
1243	GLOBAL_MATRIX.forget();
1244	}
1245	}
1246	else if (need_recalc.matrix) {
1247	GLOBAL_MATRIX.forget();
1248	}
1249
1250	error = di_calculate_matrix(aw_root, weighted_filter, bootstrap_flag, !bootstrap_flag, &aborted);
1251	if (error && aborted) {
1252	error = NULp; // clear error (otherwise no tree will be read below)
1253	break; // end of bootstrap
1254	}
1255
1256	if (!GLOBAL_MATRIX.exists()) {
1257	error = "unexpected error in di_calculate_matrix_cb (data missing)";
1258	break;
1259	}
1260
1261	DI_MATRIX *matr = GLOBAL_MATRIX.get();
1262	char *names = ARB_calloc<char>(matr->nentries+2);
1263
1264	for (size_t i=0; i<matr->nentries; i++) {
1265	names[i] = matr->entries[i]->name;
1266	}
1267	di_assert(matr->nentries == matr->matrix->size());
1268	tree = neighbourjoining(names, *matr->matrix, new SimpleRoot);
1269
1270	if (bootstrap_flag) {
1271	error = ctree->insert_tree_weighted(tree, matr->nentries, 1, false);
1272	UNCOVERED();
1273	destroy(tree); tree = NULp;
1274	// Warning: nodes stored as origin in deconstructed PARTs now point to deleted nodes.
1275
1276	loop_count++;
1277	progress->inc();
1278	if (!bootstrap_count) { // when waiting for kill
1279	int t = time(NULp);
1280	static int tlast = 0;
1281
1282	if (t>tlast) {
1283	progress->force_update();
1284	tlast = t;
1285	}
1286	}
1287	}
1288	free(names);
1289	} while (bootstrap_flag && loop_count != bootstrap_count);
1290
1291	if (!error) {
1292	if (bootstrap_flag) {
1293	tree = ctree->get_consensus_tree(error);
1294	progress->inc();
1295	if (!error) {
1296	error = GBT_is_invalid(tree);
1297	di_assert(!error);
1298	}
1299	}
1300
1301	if (!error) {
1302	GBDATA *gb_main = weighted_filter->get_gb_main();
1303	char *tree_name = aw_root->awar(AWAR_DIST_TREE_STD_NAME)->read_string();
1304	GB_begin_transaction(gb_main);
1305	error = GBT_write_tree(gb_main, tree_name, tree);
1306
1307	if (!error) {
1308	char *filter_name = AWT_get_combined_filter_name(aw_root, "dist");
1309	int transr = aw_root->awar(AWAR_DIST_CORR_TRANS)->read_int();
1310
1311	const char *comment;
1312	if (enum_trans_to_string[transr]) {
1313	comment = GBS_global_string("PRG=dnadist CORR=%s FILTER=%s PKG=ARB", enum_trans_to_string[transr], filter_name);
1314	}
1315	else {
1316	di_assert(transr == DI_TRANSFORMATION_FROM_TREE);
1317	const char *treename = aw_root->awar(AWAR_DIST_TREE_CURR_NAME)->read_char_pntr();
1318	comment = GBS_global_string("PRG=treedist (from '%s') PKG=ARB", treename);
1319	}
1320
1321	error = GBT_write_tree_remark(gb_main, tree_name, comment);
1322	free(filter_name);
1323	}
1324	error = GB_end_transaction(gb_main, error);
1325	free(tree_name);
1326	}
1327	}
1328
1329	UNCOVERED();
1330	destroy(tree);
1331
1332	// aw_status(); // remove 'abort' flag (@@@ got no equiv for arb_progress yet. really needed?)
1333
1334	if (bootstrap_flag) {
1335	if (all_names) delete all_names;
1336	GLOBAL_MATRIX.forget();
1337	}
1338	#if defined(DEBUG)
1339	else {
1340	di_assert(!all_names);
1341	}
1342	#endif // DEBUG
1343
1344	if (progress.isSet()) {
1345	progress->done();
1346	}
1347
1348	if (error) {
1349	aw_message(error);
1350	}
1351	else {
1352	need_recalc.tree = false;
1353	aw_root->awar(AWAR_TREE_REFRESH)->touch();
1354	}
1355	}
1356
1357
1358	static void di_autodetect_callback(AW_window aww, GBDATA gb_main) {
1359	GB_push_transaction(gb_main);
1360
1361	GLOBAL_MATRIX.forget();
1362
1363	AW_root *aw_root = aww->get_root();
1364	char *use = aw_root->awar(AWAR_DIST_ALIGNMENT)->read_string();
1365	long ali_len = GBT_get_alignment_len(gb_main, use);
1366	GB_ERROR error = NULp;
1367
1368	if (ali_len<=0) {
1369	GB_pop_transaction(gb_main);
1370	error = "Please select a valid alignment";
1371	GB_clear_error();
1372	}
1373	else {
1374	arb_progress progress("Analyzing data");
1375
1376	char *filter_str = aw_root->awar(AWAR_DIST_FILTER_FILTER)->read_string();
1377	// char *cancel = aw_root->awar(AWAR_DIST_CANCEL_CHARS)->read_string();
1378
1379	AliView *aliview = NULp;
1380	{
1381	AP_filter *ap_filter = NULp;
1382	long flen = strlen(filter_str);
1383
1384	if (flen == ali_len) {
1385	ap_filter = new AP_filter(filter_str, "0", ali_len);
1386	}
1387	else {
1388	if (flen) {
1389	aw_message("Warning: your filter len is not equal to the alignment len\nfilter got truncated with zeros or cutted");
1390	ap_filter = new AP_filter(filter_str, "0", ali_len);
1391	}
1392	else {
1393	ap_filter = new AP_filter(ali_len); // unfiltered
1394	}
1395	}
1396
1397	error = ap_filter->is_invalid();
1398	if (!error) {
1399	AP_weights ap_weights(ap_filter);
1400	aliview = new AliView(gb_main, *ap_filter, ap_weights, use);
1401	}
1402	delete ap_filter;
1403	}
1404
1405	if (error) {
1406	GB_pop_transaction(gb_main);
1407	}
1408	else {
1409	DI_MATRIX phm(*aliview);
1410
1411	{
1412	char *load_what = aw_root->awar(AWAR_DIST_WHICH_SPECIES)->read_string();
1413	char *sort_tree_name = aw_root->awar(AWAR_DIST_TREE_SORT_NAME)->read_string();
1414
1415	LoadWhat all_flag = (strcmp(load_what, "all") == 0) ? DI_LOAD_ALL : DI_LOAD_MARKED;
1416
1417	GLOBAL_MATRIX.forget();
1418
1419	MatrixOrder order(gb_main, sort_tree_name);
1420	error = phm.load(all_flag, order, true, NULp);
1421
1422	free(sort_tree_name);
1423	free(load_what);
1424	}
1425
1426	GB_pop_transaction(gb_main);
1427
1428	if (!error) {
1429	progress.subtitle("Search Correction");
1430
1431	string msg;
1432	DI_TRANSFORMATION detected = phm.detect_transformation(msg);
1433	aw_root->awar(AWAR_DIST_CORR_TRANS)->write_int(detected);
1434	aw_message(msg.c_str());
1435	}
1436	}
1437
1438	// free(cancel);
1439	delete aliview;
1440
1441	free(filter_str);
1442	}
1443
1444	if (error) aw_message(error);
1445
1446	free(use);
1447	}
1448
1449	__ATTR__NORETURN static void di_exit(AW_window aww, GBDATA gb_main) {
1450	if (gb_main) {
1451	AW_root *aw_root = aww->get_root();
1452	shutdown_macro_recording(aw_root);
1453	aw_root->unlink_awars_from_DB(gb_main);
1454	GB_close(gb_main);
1455	}
1456	GLOBAL_MATRIX.set_changed_cb(NULp);
1457	exit(EXIT_SUCCESS);
1458	}
1459
1460	static void di_calculate_full_matrix_cb(AW_window aww, const WeightedFilter weighted_filter) {
1461	recalculate_matrix_cb = new BoundWindowCallback(aww, makeWindowCallback(di_calculate_full_matrix_cb, weighted_filter));
1462
1463	GLOBAL_MATRIX.forget_if_not_has_type(DI_MATRIX_FULL);
1464	GB_ERROR error = di_calculate_matrix(aww->get_root(), weighted_filter, 0, true, NULp);
1465	aw_message_if(error);
1466	last_matrix_calculation_error = error;
1467	if (!error) tree_needs_recalc_cb();
1468	}
1469
1470	static void di_calculate_compressed_matrix_cb(AW_window aww, WeightedFilter weighted_filter) {
1471	recalculate_matrix_cb = new BoundWindowCallback(aww, makeWindowCallback(di_calculate_compressed_matrix_cb, weighted_filter));
1472
1473	GBDATA *gb_main = weighted_filter->get_gb_main();
1474	GB_transaction ta(gb_main);
1475
1476	AW_root *aw_root = aww->get_root();
1477	char *treename = aw_root->awar(AWAR_DIST_TREE_COMP_NAME)->read_string();
1478	GB_ERROR error = NULp;
1479	TreeNode *tree = GBT_read_tree(gb_main, treename, new SimpleRoot);
1480
1481	if (!tree) {
1482	error = GB_await_error();
1483	}
1484	else {
1485	{
1486	LocallyModify<MatrixDisplay*> skipRefresh(matrixDisplay, NULp); // skip refresh, until matrix has been compressed
1487
1488	GLOBAL_MATRIX.forget(); // always forget (as tree might have changed)
1489	error = di_calculate_matrix(aw_root, weighted_filter, 0, true, NULp);
1490	if (!error && !GLOBAL_MATRIX.exists()) {
1491	error = "Failed to calculate your matrix (bug?)";
1492	}
1493	if (!error) {
1494	error = GLOBAL_MATRIX.get()->compress(tree);
1495	}
1496	}
1497	UNCOVERED();
1498	destroy(tree);
1499
1500	// now force refresh
1501	if (matrixDisplay) {
1502	matrixDisplay->mark(MatrixDisplay::NEED_SETUP);
1503	matrixDisplay->update_display();
1504	}
1505	}
1506	free(treename);
1507	aw_message_if(error);
1508	last_matrix_calculation_error = error;
1509	if (!error) tree_needs_recalc_cb();
1510	}
1511
1512	static void di_define_sort_tree_name_cb(AW_window *aww) {
1513	AW_root *aw_root = aww->get_root();
1514	char *tree_name = aw_root->awar(AWAR_DIST_TREE_CURR_NAME)->read_string();
1515	aw_root->awar(AWAR_DIST_TREE_SORT_NAME)->write_string(tree_name);
1516	free(tree_name);
1517	}
1518	static void di_define_compression_tree_name_cb(AW_window *aww) {
1519	AW_root *aw_root = aww->get_root();
1520	char *tree_name = aw_root->awar(AWAR_DIST_TREE_CURR_NAME)->read_string();
1521	aw_root->awar(AWAR_DIST_TREE_COMP_NAME)->write_string(tree_name);
1522	free(tree_name);
1523	}
1524
1525	static void di_define_save_tree_name_cb(AW_window *aww) {
1526	AW_root *aw_root = aww->get_root();
1527	char *tree_name = aw_root->awar(AWAR_DIST_TREE_CURR_NAME)->read_string();
1528	aw_root->awar(AWAR_DIST_TREE_STD_NAME)->write_string(tree_name);
1529	free(tree_name);
1530	}
1531
1532
1533	AW_window DI_create_matrix_window(AW_root aw_root, GBDATA *gb_main) {
1534	AW_window_simple_menu *aws = new AW_window_simple_menu;
1535	aws->init(aw_root, "NEIGHBOUR JOINING", "NEIGHBOUR JOINING [ARB_DIST]");
1536	aws->load_xfig("di_ge_ma.fig");
1537	aws->button_length(10);
1538
1539	WindowCallback close_cb = makeWindowCallback(di_exit, gb_main);
1540
1541	aws->at("close");
1542	aws->callback(close_cb);
1543	aws->create_button("CLOSE", "CLOSE", "C");
1544
1545	aws->at("help");
1546	aws->callback(makeHelpCallback("dist.hlp"));
1547	aws->create_button("HELP", "HELP", "H");
1548
1549	GB_push_transaction(gb_main);
1550
1551	#if defined(DEBUG)
1552	AWT_create_debug_menu(aws);
1553	#endif // DEBUG
1554
1555	aws->create_menu("File", "F", AWM_ALL);
1556	insert_macro_menu_entry(aws, false);
1557	aws->insert_menu_topic("quit", "Quit", "Q", "quit.hlp", AWM_ALL, close_cb);
1558
1559	aws->create_menu("Properties", "P", AWM_ALL);
1560	aws->insert_menu_topic("frame_props", "Frame settings ...", "F", "props_frame.hlp", AWM_ALL, AW_preset_window);
1561	aws->sep______________();
1562	AW_insert_common_property_menu_entries(aws);
1563	aws->sep______________();
1564	aws->insert_menu_topic("save_props", "Save Properties (dist.arb)", "S", "savedef.hlp", AWM_ALL, AW_save_properties);
1565
1566	aws->insert_help_topic("ARB_DIST help", "D", "dist.hlp", AWM_ALL, makeHelpCallback("dist.hlp"));
1567
1568	// ------------------
1569	// left side
1570
1571	aws->at("which_species");
1572	aws->create_option_menu(AWAR_DIST_WHICH_SPECIES, true);
1573	aws->insert_option("all", "a", "all");
1574	aws->insert_default_option("marked", "m", "marked");
1575	aws->update_option_menu();
1576
1577	aws->at("which_alignment");
1578	awt_create_ALI_selection_list(gb_main, (AW_window )aws, AWAR_DIST_ALIGNMENT, "=");
1579
1580	// filter & weights
1581
1582	AW_awar *awar_dist_alignment = aws->get_root()->awar_string(AWAR_DIST_ALIGNMENT);
1583	WeightedFilter *weighted_filter = // do NOT free (bound to callbacks)
1584	new WeightedFilter(gb_main, aws->get_root(), AWAR_DIST_FILTER_NAME, AWAR_DIST_COLUMN_STAT_NAME, awar_dist_alignment);
1585
1586	aws->at("filter_select");
1587	aws->callback(makeCreateWindowCallback(awt_create_select_filter_win, weighted_filter->get_adfiltercbstruct()));
1588	aws->create_button("SELECT_FILTER", AWAR_DIST_FILTER_NAME);
1589
1590	aws->at("weights_select");
1591	aws->sens_mask(AWM_EXP);
1592	aws->callback(makeCreateWindowCallback(COLSTAT_create_selection_window, weighted_filter->get_column_stat()));
1593	aws->create_button("SELECT_COL_STAT", AWAR_DIST_COLUMN_STAT_NAME);
1594	aws->sens_mask(AWM_ALL);
1595
1596	aws->at("which_cancel");
1597	aws->create_input_field(AWAR_DIST_CANCEL_CHARS, 12);
1598
1599	aws->at("cancel_select");
1600	aws->callback(awt_create_select_cancel_window);
1601	aws->create_button("SELECT_CANCEL_CHARS", "Info", "C");
1602
1603	aws->at("change_matrix");
1604	aws->callback(create_dna_matrix_window);
1605	aws->create_button("EDIT_MATRIX", "Edit Matrix");
1606
1607	aws->at("enable");
1608	aws->create_toggle(AWAR_DIST_MATRIX_DNA_ENABLED);
1609
1610	aws->at("which_correction");
1611	aws->create_option_menu(AWAR_DIST_CORR_TRANS, true);
1612	aws->insert_option("none", "n", (int)DI_TRANSFORMATION_NONE);
1613	aws->insert_option("similarity", "n", (int)DI_TRANSFORMATION_SIMILARITY);
1614	aws->insert_option("jukes-cantor (dna)", "c", (int)DI_TRANSFORMATION_JUKES_CANTOR);
1615	aws->insert_option("felsenstein (dna)", "f", (int)DI_TRANSFORMATION_FELSENSTEIN);
1616	aws->insert_option("olsen (dna)", "o", (int)DI_TRANSFORMATION_OLSEN);
1617	aws->insert_option("felsenstein/voigt (exp)", "1", (int)DI_TRANSFORMATION_FELSENSTEIN_VOIGT);
1618	aws->insert_option("olsen/voigt (exp)", "2", (int)DI_TRANSFORMATION_OLSEN_VOIGT);
1619	aws->insert_option("kimura (pro)", "k", (int)DI_TRANSFORMATION_KIMURA);
1620	aws->insert_option("PAM (protein)", "c", (int)DI_TRANSFORMATION_PAM);
1621	aws->insert_option("Cat. Hall(exp)", "c", (int)DI_TRANSFORMATION_CATEGORIES_HALL);
1622	aws->insert_option("Cat. Barker(exp)", "c", (int)DI_TRANSFORMATION_CATEGORIES_BARKER);
1623	aws->insert_option("Cat.Chem (exp)", "c", (int)DI_TRANSFORMATION_CATEGORIES_CHEMICAL);
1624	aws->insert_option("from selected tree", "t", (int)DI_TRANSFORMATION_FROM_TREE);
1625	aws->insert_default_option("unknown", "u", (int)DI_TRANSFORMATION_NONE);
1626
1627	aws->update_option_menu();
1628
1629	aws->at("autodetect"); // auto
1630	aws->callback(makeWindowCallback(di_autodetect_callback, gb_main));
1631	aws->sens_mask(AWM_EXP);
1632	aws->create_button("AUTODETECT_CORRECTION", "AUTODETECT", "A");
1633	aws->sens_mask(AWM_ALL);
1634
1635	// -------------------
1636	// right side
1637
1638
1639	aws->at("mark_distance");
1640	aws->callback(makeWindowCallback(di_mark_by_distance, weighted_filter));
1641	aws->create_autosize_button("MARK_BY_DIST", "Mark all species");
1642
1643	aws->at("mark_lower");
1644	aws->create_input_field(AWAR_DIST_MIN_DIST, 5);
1645
1646	aws->at("mark_upper");
1647	aws->create_input_field(AWAR_DIST_MAX_DIST, 5);
1648
1649	// -----------------
1650
1651	// tree selection
1652
1653	aws->at("tree_list");
1654	awt_create_TREE_selection_list(gb_main, aws, AWAR_DIST_TREE_CURR_NAME, true);
1655
1656	aws->at("detect_clusters");
1657	aws->callback(makeCreateWindowCallback(DI_create_cluster_detection_window, weighted_filter));
1658	aws->create_autosize_button("DETECT_CLUSTERS", "Detect homogenous clusters in tree", "D");
1659
1660	// matrix calculation
1661
1662	aws->button_length(18);
1663
1664	aws->at("calculate");
1665	aws->callback(makeWindowCallback(di_calculate_full_matrix_cb, weighted_filter));
1666	aws->create_button("CALC_FULL_MATRIX", "Calculate\nFull Matrix", "F");
1667
1668	aws->at("compress");
1669	aws->callback(makeWindowCallback(di_calculate_compressed_matrix_cb, weighted_filter));
1670	aws->create_button("CALC_COMPRESSED_MATRIX", "Calculate\nCompressed Matrix", "C");
1671
1672	recalculate_matrix_cb = new BoundWindowCallback(aws, makeWindowCallback(di_calculate_full_matrix_cb, weighted_filter));
1673
1674	aws->button_length(13);
1675
1676	{
1677	static save_matrix_params sparams;
1678
1679	sparams.awar_base = AWAR_DIST_SAVE_MATRIX_BASE;
1680	sparams.weighted_filter = weighted_filter;
1681
1682	aws->at("save_matrix");
1683	aws->callback(makeCreateWindowCallback(DI_create_save_matrix_window, &sparams));
1684	aws->create_button("SAVE_MATRIX", "Save matrix", "M");
1685
1686	aws->at("view_matrix");
1687	aws->callback(makeWindowCallback(di_view_matrix_cb, &sparams));
1688	aws->create_button("VIEW_MATRIX", "View matrix", "V");
1689	}
1690
1691	aws->button_length(22);
1692	aws->at("use_compr_tree");
1693	aws->callback(di_define_compression_tree_name_cb);
1694	aws->create_button("USE_COMPRESSION_TREE", "Use to compress", "");
1695	aws->at("use_sort_tree");
1696	aws->callback(di_define_sort_tree_name_cb);
1697	aws->create_button("USE_SORT_TREE", "Use to sort", "");
1698
1699	aws->at("compr_tree_name"); aws->create_input_field(AWAR_DIST_TREE_COMP_NAME, 12);
1700	aws->at("sort_tree_name"); aws->create_input_field(AWAR_DIST_TREE_SORT_NAME, 12);
1701
1702	// tree calculation
1703
1704	aws->button_length(18);
1705
1706	aws->at("t_calculate");
1707	aws->callback(makeWindowCallback(di_calculate_tree_cb, weighted_filter, false));
1708	aws->create_button("CALC_TREE", "Calculate \ntree", "C");
1709
1710	aws->at("bootstrap");
1711	aws->callback(makeWindowCallback(di_calculate_tree_cb, weighted_filter, true));
1712	aws->create_button("CALC_BOOTSTRAP_TREE", "Calculate \nbootstrap tree");
1713
1714	recalculate_tree_cb = new BoundWindowCallback(aws, makeWindowCallback(di_calculate_tree_cb, weighted_filter, false));
1715
1716	aws->button_length(22);
1717	aws->at("use_existing");
1718	aws->callback(di_define_save_tree_name_cb);
1719	aws->create_button("USE_NAME", "Use as new tree name", "");
1720
1721	aws->at("calc_tree_name");
1722	aws->create_input_field(AWAR_DIST_TREE_STD_NAME, 12);
1723
1724	aws->at("bcount");
1725	aws->create_input_field(AWAR_DIST_BOOTSTRAP_COUNT, 7);
1726
1727	{
1728	aws->sens_mask(AWM_EXP);
1729
1730	aws->at("auto_calc_tree");
1731	aws->label("Auto calculate tree");
1732	aws->create_toggle(AWAR_DIST_MATRIX_AUTO_CALC_TREE);
1733
1734	aws->at("auto_recalc");
1735	aws->label("Auto recalculate");
1736	aws->create_toggle(AWAR_DIST_MATRIX_AUTO_RECALC);
1737
1738	aws->sens_mask(AWM_ALL);
1739	}
1740
1741	bool disable_autocalc = !ARB_in_expert_mode(aw_root);
1742	if (disable_autocalc) {
1743	aw_root->awar(AWAR_DIST_MATRIX_AUTO_RECALC)->write_int(0);
1744	aw_root->awar(AWAR_DIST_MATRIX_AUTO_CALC_TREE)->write_int(0);
1745	}
1746
1747	GB_pop_transaction(gb_main);
1748	return aws;
1749	}
1750
1751	// --------------------------------------------------------------------------------
1752
1753	#ifdef UNIT_TESTS
1754	#include <arb_diff.h>
1755	#include <arb_file.h>
1756
1757	#ifndef TEST_UNIT_H
1758	#include <test_unit.h>
1759	#endif
1760
1761	class DIST_testenv : virtual Noncopyable {
1762	GB_shell shell;
1763	GBDATA *gb_main;
1764	AliView *ali_view;
1765
1766	public:
1767	DIST_testenv(const char dbname, const char aliName)
1768	: ali_view(NULp)
1769	{
1770	gb_main = GB_open(dbname, "r");
1771	TEST_REJECT_NULL(gb_main);
1772
1773	GB_transaction ta(gb_main);
1774	size_t aliLength = GBT_get_alignment_len(gb_main, aliName);
1775
1776	AP_filter filter(aliLength);
1777	if (!filter.is_invalid()) {
1778	AP_weights weights(&filter);
1779	ali_view = new AliView(gb_main, filter, weights, aliName);
1780	}
1781	}
1782	~DIST_testenv() {
1783	delete ali_view;
1784	GB_close(gb_main);
1785	}
1786
1787	const AliView& aliview() const { return *ali_view; }
1788	GBDATA *gbmain() const { return gb_main; }
1789	};
1790
1791	void TEST_matrix() {
1792	for (int iat = GB_AT_RNA; iat<=GB_AT_AA; ++iat) {
1793	GB_alignment_type at = GB_alignment_type(iat);
1794	// ---------------
1795	// setup
1796	// GB_AT_RNA GB_AT_DNA GB_AT_AA
1797	const char *db_name[]= { NULp, "TEST_trees.arb", "TEST_realign.arb", "TEST_realign.arb", };
1798	const char *ali_name[]= { NULp, "ali_5s", "ali_dna", "ali_pro", };
1799
1800	TEST_ANNOTATE(GBS_global_string("ali_name=%s", ali_name[at]));
1801	DIST_testenv env(db_name[at], ali_name[at]);
1802
1803	DI_MATRIX matrix(env.aliview());
1804	MatrixOrder order(env.gbmain(), "tree_abc"); // no such tree!
1805	TEST_EXPECT_NO_ERROR(matrix.load(DI_LOAD_MARKED, order, true, NULp));
1806
1807	// -------------------------------
1808	// detect_transformation
1809	DI_TRANSFORMATION detected_trans;
1810	{
1811	string msg;
1812
1813	detected_trans = matrix.detect_transformation(msg);
1814	DI_TRANSFORMATION expected = DI_TRANSFORMATION_NONE_DETECTED;
1815	switch (at) {
1816	case GB_AT_RNA: expected = DI_TRANSFORMATION_NONE; break;
1817	case GB_AT_DNA: expected = DI_TRANSFORMATION_JUKES_CANTOR; break;
1818	case GB_AT_AA: expected = DI_TRANSFORMATION_PAM; break;
1819	case GB_AT_UNKNOWN: di_assert(0); break;
1820	}
1821	TEST_EXPECT_EQUAL(detected_trans, expected);
1822	}
1823
1824	// ------------------------------
1825	// calculate the matrix
1826
1827	// @@@ does not test user-defined transformation-matrix!
1828	if (at == GB_AT_AA) {
1829	matrix.calculate_pro(detected_trans, NULp);
1830	}
1831	else {
1832	if (at == GB_AT_RNA) detected_trans = DI_TRANSFORMATION_FELSENSTEIN; // force calculate_overall_freqs
1833	matrix.calculate("", detected_trans, NULp, NULp);
1834	}
1835
1836	// -----------------------------------
1837	// save in available formats
1838
1839	for (DI_SAVE_TYPE saveType = DI_SAVE_PHYLIP_COMP; saveType<=DI_SAVE_TABBED; saveType = DI_SAVE_TYPE(saveType+1)) {
1840	const char *savename = "distance/matrix.out";
1841	matrix.save(savename, saveType);
1842
1843	const char *suffixAT[] = { NULp, "rna", "dna", "pro" };
1844	const char *suffixST[] = { "phylipComp", "readable", "tabbed" };
1845	char *expected = GBS_global_string_copy("distance/matrix.%s.%s.expected", suffixAT[at], suffixST[saveType]);
1846
1847	// #define TEST_AUTO_UPDATE // uncomment to auto-update expected matrices
1848
1849	#if defined(TEST_AUTO_UPDATE)
1850	TEST_COPY_FILE(savename, expected);
1851	#else
1852	TEST_EXPECT_TEXTFILES_EQUAL(savename, expected);
1853	#endif // TEST_AUTO_UPDATE
1854	TEST_EXPECT_ZERO_OR_SHOW_ERRNO(GB_unlink(savename));
1855
1856	free(expected);
1857	}
1858	}
1859	}
1860
1861	#endif // UNIT_TESTS
1862
1863	// --------------------------------------------------------------------------------
1864

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

Download in other formats: