Context Navigation

source: tags/ms_ra2q56/DIST/DI_matr.cxx

Visit:

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

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

Download in other formats: