1/*---------------------------------------------------------*/
2/* Turn the keys of an AVL into a string for labeling
3   chart rows */
4
5
6
7function avlToLabels (_gb_anAVL,_gb_prefix,_gb_delim)
8{
9	_gb_resString = "";
10	_gb_keys	  = Rows (_gb_anAVL);
11	_gb_count	  = Columns (_gb_keys);
12	_gb_resString * 128;
13	_gb_resString * _gb_prefix;
14	if (Abs(_gb_prefix))
15	{
16		_gb_resString * _gb_delim;
17	}
18	if (_gb_count)
19	{
20		_gb_resString * _gb_keys[0];
21	}
22	for (_gb_idx = 1; _gb_idx < _gb_count; _gb_idx = _gb_idx + 1)
23	{
24		_gb_resString * (_gb_delim+_gb_keys[_gb_idx]);
25	}
26	_gb_resString * 0;
27	return _gb_resString;
28}
29
30/*---------------------------------------------------------*/
31/* Turn the keys of an AVL into a numerical column matrix  */
32
33function avlKeysToMatrix (_gb_anAVL)
34{
35	_gb_keys	  = Rows (_gb_anAVL);
36	_gb_count	  = Columns (_gb_keys);
37	_gb_resMatrix = {_gb_count,1};
38
39	for (_gb_idx = 0; _gb_idx < _gb_count; _gb_idx = _gb_idx + 1)
40	{
41		_gb_resMatrix[_gb_idx] = 0+_gb_keys[_gb_idx];
42	}
43	return _gb_resMatrix;
44}
45
46/*---------------------------------------------------------*/
47/* Assuming that the AVL is 0..N indexed, produce a
48string with AVL entries separated by _gb_delim */
49
50function avlToString (_gb_anAVL,_gb_delim)
51{
52	_gb_count	  = Abs (_gb_anAVL);
53	_gb_resString = "";
54	_gb_resString * 128;
55	if (_gb_count)
56	{
57		_gb_resString * (""+_gb_anAVL[0]);
58	}
59	for (_gb_idx = 1; _gb_idx < _gb_count; _gb_idx = _gb_idx + 1)
60	{
61		_gb_resString * (_gb_delim+_gb_anAVL[_gb_idx]);
62	}
63	_gb_resString * 0;
64	return _gb_resString;
65}
66
67
68/*---------------------------------------------------------*/
69/* Stratify avl by values; store indices for each unique value */
70
71function stratifyAVLByValuesAux (key,value)
72{
73	if (Abs (_gb_resAVL[value]) == 0)
74	{
75		_gb_resAVL[value] = {};
76	}
77
78	_gb_resAVL[value] + key;
79
80	return 0;
81
82}
83
84function stratifyAVLByValues (_gb_anAVL)
85{
86	_gb_count	  = Abs (_gb_anAVL);
87	_gb_resAVL	  = {};
88	_gb_anAVL["stratifyAVLByValuesAux"][""];
89	return _gb_resAVL;
90}
91
92/*---------------------------------------------------------*/
93/* Stratify a matrix by values; store indices for each unique value */
94
95function stratifyMatrixByValues (_gb_aMatrix)
96{
97	_gb_count	  = Rows (_gb_aMatrix)*Columns(_gb_aMatrix);
98	_gb_resAVL	  = {};
99	for (_gb_idx = 0; _gb_idx < _gb_count; _gb_idx += 1)
100	{
101		_gb_key = _gb_aMatrix[_gb_idx];
102		if (Abs (_gb_resAVL[_gb_key]) == 0)
103		{
104			_gb_resAVL[_gb_key] = {};
105		}
106		(_gb_resAVL[_gb_key]) + _gb_idx;
107	}
108	return _gb_resAVL;
109}
110
111
112/*---------------------------------------------------------*/
113/* Assuming that the AVL is 0..N indexed, produce a
114row matrix with AVL entries, using _gb_map to map the values
115and _gb_stride to do the conversion */
116
117function avlToRow (_gb_anAVL, _gb_map, _gb_stride)
118{
119	_gb_count	  = Abs (_gb_anAVL);
120	_gb_matrix	  = {1,_gb_count*_gb_stride};
121
122	if (_gb_stride>1)
123	{
124		for (_gb_idx = 0; _gb_idx < _gb_count; _gb_idx = _gb_idx + 1)
125		{
126			for (_gb_idx2 = 0; _gb_idx2 < _gb_stride; _gb_idx2 = _gb_idx2 + 1)
127			{
128				_gb_matrix [_gb_idx*_gb_stride+_gb_idx2] = _gb_map[_gb_stride*_gb_anAVL[_gb_idx]+_gb_idx2];
129			}
130		}
131	}
132	else
133	{
134		for (_gb_idx = 0; _gb_idx < _gb_count; _gb_idx = _gb_idx + 1)
135		{
136			_gb_matrix [_gb_idx] = _gb_map[_gb_anAVL[_gb_idx]];
137		}
138	}
139	return _gb_matrix;
140}
141
142/*---------------------------------------------------------*/
143/* Split a file path into DIRECTORY, FILENAME and EXTENSION  */
144
145function splitFilePath (_filePath)
146{
147	_splitPath = {};
148	_split     = _filePath $ ("[^\\"+DIRECTORY_SEPARATOR+"]+$");
149	if (_split[0] == 0 && _split[1] == Abs (_filePath)-1) /* no path, all file name */
150	{
151		_splitPath ["DIRECTORY"] = "";
152	}
153	else
154	{
155		_splitPath ["DIRECTORY"] = _filePath[0][_split[0]-1];
156		_filePath = _filePath[_split[0]][Abs(_filePath)];
157	}
158
159	_split     = _filePath || "\\.";
160	if (_split[0] < 0) /* no extension */
161	{
162		_splitPath ["EXTENSION"] = "";
163		_splitPath ["FILENAME"]  = _filePath;
164 	}
165	else
166	{
167		_splitPath ["EXTENSION"] = _filePath[_split[Rows(_split)-1]+1][Abs(_filePath)-1];
168		_splitPath ["FILENAME"]  = _filePath[0][_split[Rows(_split)-1]-1];
169	}
170	return _splitPath;
171}
172
173/*---------------------------------------------------------*/
174/* fix global variables in a LF at their current values */
175
176function fixGlobalParameters (_lfName) {
177	GetString (_lfInfo,^_lfName,-1);
178	_lfInfo = _lfInfo["Global Independent"];
179	for (_gb_idx = 0; _gb_idx < Columns (_lfInfo); _gb_idx += 1) {
180		ExecuteCommands (_lfInfo[_gb_idx] + ":=" + _lfInfo[_gb_idx] + "__;");
181	}
182	return 0;
183}
184
185/*---------------------------------------------------------*/
186/* unconstrain global variables in a LF at their current values */
187
188function unconstrainGlobalParameters (_lfName) {
189	GetString (_lfInfo,^_lfName,-1);
190	_lfInfo = _lfInfo["Global Constrained"];
191	for (_gb_idx = 0; _gb_idx < Columns (_lfInfo); _gb_idx += 1) {
192		ExecuteCommands (_lfInfo[_gb_idx] + "=" + _lfInfo[_gb_idx]);
193	}
194	return 0;
195}
196
197/*---------------------------------------------------------*/
198/* prompt for global variabless in a LF and fix their values */
199
200function promptForGlobalParameters (_lfName)
201{
202	ExecuteCommands ("GetString (_lfInfo," + _lfName + ",-1);");
203	_lfInfo = _lfInfo["Global Independent"];
204	for (_gb_idx = 0; _gb_idx < Columns (_lfInfo); _gb_idx = _gb_idx + 1)
205	{
206		fprintf (stdout, "\nEnter a value for ", _lfInfo[_gb_idx], ":");
207		fscanf  (stdin, "Number", _tval);
208		ExecuteCommands (_lfInfo[_gb_idx] + ":=" + _tval + ";");
209	}
210	return 0;
211}
212
213/*---------------------------------------------------------*/
214/* echo global parameters */
215
216function echoGlobalParameters (_lfName)
217{
218	ExecuteCommands ("GetString (_lfInfo," + _lfName + ",-1);");
219	_lfInfo = _lfInfo["Global Independent"];
220	for (_gb_idx = 0; _gb_idx < Columns (_lfInfo); _gb_idx = _gb_idx + 1)
221	{
222		ExecuteCommands ("_tval = "+_lfInfo[_gb_idx]);
223		fprintf (stdout, _lfInfo[_gb_idx], " : ", Format (_tval, 12, 4), "\n");
224	}
225	return Columns (_lfInfo);
226}
227
228
229/*---------------------------------------------------------*/
230/* take a snapshot of global parameters */
231
232function stashGlobalParameters (_lfName)
233{
234	ExecuteCommands ("GetString (_lfInfo," + _lfName + ",-1);");
235	_lfInfo = _lfInfo["Global Independent"];
236	_paramStash = {};
237	for (_gb_idx = 0; _gb_idx < Columns (_lfInfo); _gb_idx = _gb_idx + 1)
238	{
239		ExecuteCommands ("_paramStash[\""+_lfInfo[_gb_idx]+"\"] =" + _lfInfo[_gb_idx] + ";");
240	}
241	return _paramStash;
242}
243
244/*---------------------------------------------------------*/
245/* define a global parameter if not already defined */
246
247function defineIfNeeded (_parName, _parValue)
248{
249	ExecuteCommands("GetInformation (_gb_idx, \"^`_parName`$\");");
250	if (Rows (_gb_idx) == 0)
251	{
252		ExecuteCommands ("global `_parName`="+_parValue+";");
253		return 0;
254	}
255	return 1;
256}
257
258/*---------------------------------------------------------*/
259/* export a list of variables into the string of the form
260
261   varID = varValue;
262
263*/
264
265function exportVarList (_varList)
266{
267	_exportString = ""; _exportString * 256;
268
269	for (_idx = 0; _idx < Columns (_varList) * Rows (_varList); _idx += 1)
270	{
271		_exportString * (_varList [_idx] + " = " + Eval (_varList[_idx]) + ";\n");
272	}
273
274	_exportString * 0;
275	return _exportString;
276}
277
278/*---------------------------------------------------------*/
279/* restore values of global parameters */
280
281function restoreGlobalParameters (_paramStash)
282{
283	_stashKeys = Rows(_paramStash);
284	for (_gb_idx = 0; _gb_idx < Abs (_paramStash); _gb_idx = _gb_idx + 1)
285	{
286		ExecuteCommands (_stashKeys[_gb_idx] + "=" + _paramStash[_stashKeys[_gb_idx]] + ";");
287	}
288	return 0;
289}
290
291/*---------------------------------------------------------*/
292/* take a string row/columqn matrix and turn it into an AVL of
293   the form avl["matrix entry"] = 1 for each matrix entry */
294
295function stringMatrixToAVL (_theList&)
296{
297	_gb_dim = Rows(_theList)*Columns(_theList);
298	_gb_ret = {};
299	for (_gb_idx = 0; _gb_idx < _gb_dim; _gb_idx = _gb_idx + 1)
300	{
301		_gb_ret [_theList[_gb_idx]] = _gb_idx+1;
302	}
303	return _gb_ret;
304}
305
306/*---------------------------------------------------------*/
307/* take an avl indexed by 0..N-1 and convert to a row matrix */
308
309function avlToMatrix (_theList&)
310{
311	_gb_dim = Abs(_theList);
312	_gb_ret = {_gb_dim,1};
313	for (_gb_idx = 0; _gb_idx < _gb_dim; _gb_idx = _gb_idx + 1)
314	{
315		_gb_ret [_gb_idx] = _theList[_gb_idx];
316	}
317	return _gb_ret;
318}
319
320/*---------------------------------------------------------*/
321/* take an avl indexed by 0..N-1 and convert to a choice list matrix */
322
323function avlToChoiceMatrix (_theList&)
324{
325	_gb_dim = Abs(_theList);
326	_gb_ret = {_gb_dim,2};
327	for (_gb_idx = 0; _gb_idx < _gb_dim; _gb_idx = _gb_idx + 1)
328	{
329		_gb_ret [_gb_idx][0] = _theList[_gb_idx];
330		_gb_ret [_gb_idx][1] = _theList[_gb_idx];
331	}
332	return _gb_ret;
333}
334
335/*---------------------------------------------------------*/
336/* swap the roles of integer keys (+1) and string values in an avl */
337
338function swapKeysAndValues (_theList&)
339{
340	_gb_ret 	= 	{};
341	_gb_dim 	= 	Abs(_theList);
342	_gb_keys	= 	Rows (_theList);
343
344	for (_gb_idx = 0; _gb_idx < _gb_dim; _gb_idx = _gb_idx + 1)
345	{
346		_gb_ret [_theList[_gb_keys[_gb_idx]]] = 1+_gb_keys[_gb_idx];
347	}
348	return _gb_ret;
349}
350
351
352/*---------------------------------------------------------*/
353/* report a string version of an a/b ratio, handling the cases
354   of a and/or b = 0 */
355
356function _standardizeRatio (_num, _denom)
357{
358	if (_denom != 0)
359	{
360		_ratio = _num/_denom;
361		if (_ratio < 10000)
362		{
363			return Format (_ratio,10,4);
364		}
365	}
366	else
367	{
368		if (_num == 0)
369		{
370			return " Undefined";
371		}
372	}
373	return "  Infinite";
374}
375
376/*---------------------------------------------------------*/
377/* copy branch lengths */
378
379function _copyBranchLengths (_treeID1, _treeID2, _op, _suffix)
380{
381	ExecuteCommands ("_gb_dim=BranchName("+_treeID2+",-1);");
382	_gb_ret = "";
383	_gb_ret * 128;
384
385	for (_gb_idx = 0; _gb_idx < Columns(_gb_dim)-1; _gb_idx = _gb_idx + 1)
386	{
387		_gb_idx2 = _treeID2 + "." + _gb_dim[_gb_idx] + "." + _suffix;
388		ExecuteCommands ("_gb_idx2="+_gb_idx2);
389		_gb_ret * (_treeID1 + "." + _gb_dim[_gb_idx] + "." +_suffix + ":=" + _op + _gb_idx2 + ";\n");
390	}
391	_gb_ret * 0;
392	return _gb_ret;
393}
394
395
396/*---------------------------------------------*/
397/* convert a number into a 3 letter string
398   for initializing stdin lists */
399/*---------------------------------------------*/
400
401function _mapNumberToString (n)
402{
403	if (n>=100)
404	{
405		return "" + n;
406	}
407	if (n>=10)
408	{
409		return "0" + n;
410	}
411	return "00" + n;
412}
413
414/*---------------------------------------------
415 given two integer vectors of the same length,
416cross-tabulate the values of each vector (vec1
417gives the row index, vec2 - the column index)
418---------------------------------------------*/
419
420function table (vec1, vec2)
421{
422	if (Rows(vec1) == Rows(vec2) && Columns (vec1) == 1 && Columns (vec2) == 1)
423	{
424		maxV = 0;
425		rd = Rows(vec1);
426		for (_r = 0; _r < rd; _r = _r+1)
427		{
428			maxV = Max (maxV, Max(vec1[_r],vec2[_r]));
429		}
430		tableOut = {maxV+1, maxV+1};
431		for (_r = 0; _r < rd; _r = _r+1)
432		{
433			tableOut[vec1[_r]][vec2[_r]] = tableOut[vec1[_r]][vec2[_r]] + 1;
434		}
435		return tableOut;
436	}
437	else
438	{
439		return 0;
440	}
441}
442
443/*---------------------------------------------
444 given two integer vectors of the same length,
445cross-tabulate the values of each vector (vec1
446gives the row index, vec2 - the column index)
447---------------------------------------------*/
448
449function vector_max (vec1)
450{
451	rd   = Rows(vec1) * Columns (vec1);
452	maxV = -1e100;
453	maxI = 0;
454	for (_r = 0; _r < rd; _r = _r+1)
455	{
456		if (vec1[_r] > maxV)
457		{
458			maxV = vec1[_r];
459			maxI = _r;
460		}
461	}
462	return {{maxV__,maxI__}};
463}
464
465
466/*---------------------------------------------
467 the analog of Python's string.join (list)
468---------------------------------------------*/
469
470function string_join (sep, list)
471{
472	_result = ""; _result * 128;
473
474	if (Type (list) == "Matrix")
475	{
476		_dim = Rows(list)*Columns(list);
477	}
478	else
479	{
480		_dim = Abs (list);
481	}
482
483	if (_dim)
484	{
485		_result * ("" + list[0]);
486		for (_r = 1; _r < _dim; _r = _r + 1)
487		{
488			_result * (sep + list[_r]);
489
490		}
491	}
492
493	_result * 0;
494	return _result;
495}
496
497/*---------------------------------------------
498 prompt for a value in a given range, given
499 a default value.
500---------------------------------------------*/
501
502function prompt_for_a_value (prompt,default,lowerB,upperB,isInteger)
503{
504	value = lowerB-1;
505	while (value < lowerB || value > upperB)
506	{
507		fprintf (stdout, "<<", prompt, " (permissible range = [", lowerB, ",", upperB, "], default value = ", default);
508		if (isInteger)
509		{
510			fprintf (stdout, ", integer");
511		}
512		fprintf (stdout, ")>>");
513		fscanf  (stdin, "String", strVal);
514		if (Abs(strVal) == 0)
515		{
516			value = 0+default;
517			break;
518		}
519		value = 0+strVal;
520	}
521	if (isInteger)
522		return value$1;
523
524	return value;
525}
526
527/*---------------------------------------------
528take an AVL of the form ["string"] = number
529and print it as:
530
531key[_sepChar]+: number (%)
532
533---------------------------------------------*/
534
535function _printAnAVL (_theList, _sepChar)
536{
537	_gb_keys 		= _sortStrings(Rows (_theList));
538	_printAnAVLInt (_theList, _gb_keys, _sepChar, 0);
539
540	return 0;
541}
542
543/*---------------------------------------------
544take an AVL of the form ["string"] = number
545and print it as:
546
547key[_sepChar]+: number (%)
548
549add a "Total" row
550
551---------------------------------------------*/
552
553function _printAnAVLTotal (_theList, _sepChar)
554{
555	_gb_keys 		= _sortStrings(Rows (_theList));
556	_printAnAVLInt (_theList, _gb_keys, _sepChar, 1);
557
558	return 0;
559}
560
561/*---------------------------------------------
562take an AVL of the form [number] = number
563and print it as:
564
565key[_sepChar]+: number (%)
566
567---------------------------------------------*/
568
569function _printAnAVLNumeric (_theList, _sepChar)
570{
571	_gb_dim   		= Abs(_theList);
572	num_keys		= avlKeysToMatrix (_theList)%0;
573	_gb_keys 		= {_gb_dim,1};
574	for (_gb_idx = 0; _gb_idx < _gb_dim; _gb_idx = _gb_idx + 1)
575	{
576		_gb_keys[_gb_idx] = ""+num_keys[_gb_idx];
577	}
578
579	_printAnAVLInt (_theList, _gb_keys, _sepChar, 0);
580
581	return 0;
582}
583
584function _printAnAVLNumericTotal (_theList, _sepChar)
585{
586	_gb_dim   		= Abs(_theList);
587	num_keys		= avlKeysToMatrix (_theList)%0;
588	_gb_keys 		= {_gb_dim,1};
589	for (_gb_idx = 0; _gb_idx < _gb_dim; _gb_idx = _gb_idx + 1)
590	{
591		_gb_keys[_gb_idx] = ""+num_keys[_gb_idx];
592	}
593
594	_printAnAVLInt (_theList, _gb_keys, _sepChar, 1);
595
596	return 0;
597}
598
599/*---------------------------------------------*/
600
601function _printAnAVLInt (_theList, _gb_keys, _sepChar, _doTotal)
602{
603	_gb_dim   		= Abs(_theList);
604	_gb_total 		= 0;
605	_gb_max_key_len = 0;
606	if (_doTotal)
607	{
608		gb_max_key_len = 5;
609	}
610
611	for (_gb_idx = 0; _gb_idx < _gb_dim; _gb_idx = _gb_idx + 1)
612	{
613		_gb_key 		= _gb_keys[_gb_idx];
614		_gb_max_key_len = Max (_gb_max_key_len, Abs(_gb_key));
615		_gb_total 		= _gb_total + _theList[_gb_key];
616	}
617
618	fprintf (stdout, "\n");
619	for (_gb_idx = 0; _gb_idx < _gb_dim; _gb_idx = _gb_idx + 1)
620	{
621		_gb_key 		= _gb_keys[_gb_idx];
622		fprintf (stdout, _gb_key);
623		for (_gb_idx2 = Abs(_gb_key); _gb_idx2 < _gb_max_key_len; _gb_idx2 = _gb_idx2 + 1)
624		{
625			fprintf (stdout, _sepChar);
626		}
627		fprintf (stdout, ":", Format (_theList[_gb_key],8,0), " (", Format (100*_theList[_gb_key]/_gb_total,5,2), "%)\n");
628	}
629
630	if (_doTotal)
631	{
632		fprintf (stdout, "Total");
633		for (_gb_idx2 = 5; _gb_idx2 < _gb_max_key_len; _gb_idx2 = _gb_idx2 + 1)
634		{
635			fprintf (stdout, _sepChar);
636		}
637		fprintf (stdout, ":", Format (_gb_total,8,0),"\n");
638	}
639
640	return 0;
641}
642
643/*---------------------------------------------
644sort a matrix of strings; return a
645column vector
646---------------------------------------------*/
647function _sortStringsAux (theKey, theValue)
648{
649	for (_gb_idx2 = 0; _gb_idx2 < theValue; _gb_idx2=_gb_idx2+1)
650	{
651		_gb_sortedStrings [_gb_idx] = theKey;
652		_gb_idx = _gb_idx + 1;
653	}
654	return 0;
655}
656
657function _sortStrings (_theList)
658{
659	_gb_dim = Rows (_theList)*Columns (_theList);
660	_toSort = {};
661	for (_gb_idx = 0; _gb_idx < _gb_dim; _gb_idx = _gb_idx + 1)
662	{
663		_toSort[_theList[_gb_idx]] = _toSort[_theList[_gb_idx]]+1;
664	}
665	_gb_sortedStrings = {_gb_dim,1};
666	_gb_idx = 0;
667	_toSort["_sortStringsAux"][""];
668	return _gb_sortedStrings;
669}
670
671/*---------------------------------------------
672construct the frequency vector and a LF mixing componenent for
673a general discrete distribution on numberOfRates rates
674---------------------------------------------*/
675
676function generate_gdd_freqs (numberOfRates, freqs&, lfMixing&, probPrefix, incrementFlag)
677{
678	freqs    	 = {numberOfRates,1};
679	lfMixing	 = ""; lfMixing * 128; lfMixing * "Log(";
680
681	if (numberOfRates == 1)
682	{
683		freqs[0] = "1";
684	}
685	else
686	{
687		if (incrementFlag)
688		{
689			mi = numberOfRates-1;
690			ExecuteCommands ("global "+probPrefix+"_"+mi+":<1;"+probPrefix+"_"+mi+"=1/2;");
691			ExecuteCommands ("global "+probPrefix+"_"+mi+":>0;");
692			for (mi=1; mi<numberOfRates-1; mi=mi+1)
693			{
694				ExecuteCommands (""+probPrefix+"_"+mi+" = "+probPrefix+"_"+mi+"*(1-1/numberOfRates);");
695			}
696
697		}
698		else
699		{
700			for (mi=1; mi<numberOfRates; mi=mi+1)
701			{
702				ExecuteCommands ("global "+probPrefix+"_"+mi+":<1;"+probPrefix+"_"+mi+" = 1/" + (numberOfRates-mi+1));
703				ExecuteCommands ("global "+probPrefix+"_"+mi+":>0;");
704			}
705		}
706
707		freqs[0] 	 = ""+probPrefix+"_1";
708		for (mi=1; mi<numberOfRates-1; mi=mi+1)
709		{
710			freqs[mi] = "";
711			for (mi2=1;mi2<=mi;mi2=mi2+1)
712			{
713				freqs[mi] = freqs[mi]+"(1-"+probPrefix+"_"+mi2+")";
714			}
715			freqs[mi] = freqs[mi]+""+probPrefix+"_"+(mi+1);
716		}
717
718		freqs[mi] = "";
719		for (mi2=1;mi2<mi;mi2=mi2+1)
720		{
721			freqs[mi] = freqs[mi]+"(1-"+probPrefix+"_"+mi2+")";
722		}
723		freqs[mi] = freqs[mi]+"(1-"+probPrefix+"_"+mi+")";
724	}
725
726	lfMixing * ("SITE_LIKELIHOOD[0]*"+freqs[0]);
727	for (mi = 1; mi < numberOfRates; mi=mi+1)
728	{
729		lfMixing * ("+SITE_LIKELIHOOD[" + mi + "]*" + freqs[mi]);
730	}
731	lfMixing * ")";
732	lfMixing * 0;
733	return 0;
734}
735
736/*---------------------------------------------
737reverse complement a nucleotide string
738---------------------------------------------*/
739
740_nucleotide_rc = {};
741_nucleotide_rc["A"] = "T";
742_nucleotide_rc["C"] = "G";
743_nucleotide_rc["G"] = "C";
744_nucleotide_rc["T"] = "A";
745_nucleotide_rc["M"] = "K";
746_nucleotide_rc["R"] = "Y";
747_nucleotide_rc["W"] = "W";
748_nucleotide_rc["S"] = "S";
749_nucleotide_rc["Y"] = "R";
750_nucleotide_rc["K"] = "M";
751_nucleotide_rc["B"] = "V";  /* not A */
752_nucleotide_rc["D"] = "H";  /* not C */
753_nucleotide_rc["H"] = "D";  /* not G */
754_nucleotide_rc["V"] = "B";  /* not T */
755_nucleotide_rc["N"] = "N";
756
757function nucleotideReverseComplement (seqIn)
758{
759	_seqOut = "";_seqOut*128;
760	_seqL   = Abs(seqIn);
761	for (_r = _seqL-1; _r >=0 ; _r = _r-1)
762	{
763		_seqOut *_nucleotide_rc[seqIn[_r]];
764	}
765	_seqOut*0;
766	return _seqOut;
767}
768
769/*---------------------------------------------------------------------*/
770
771function RankMatrix (matrix)
772/* take a sorted row matrix and return a rank row matrix averaging ranks for tied values */
773{
774	lastValue				   			 = matrix[0];
775	lastIndex				   			 = 0;
776	matrix							 [0] = 0;
777
778	sampleCount = Rows (matrix);
779
780	for (_i = 1; _i < sampleCount; _i = _i+1)
781	{
782		if (lastValue != matrix[_i])
783		{
784			meanIndex = _i - lastIndex;
785			lastValue = matrix[_i];
786			if (meanIndex > 1)
787			{
788				meanIndex = (lastIndex + _i - 1) * meanIndex / (2 * meanIndex);
789				for (_j = lastIndex; _j < _i; _j = _j + 1)
790				{
791					matrix[_j] = meanIndex;
792				}
793			}
794			matrix[_i] = _i;
795			lastIndex = _i;
796		}
797	}
798
799	meanIndex = _i - lastIndex;
800	if (meanIndex > 1)
801	{
802		meanIndex = (lastIndex + _i - 1) * meanIndex / (2 * meanIndex);
803		for (_j = lastIndex; _j < _i; _j = _j + 1)
804		{
805			matrix[_j] = meanIndex;
806		}
807	}
808	else
809	{
810		matrix[_i-1] = _i - 1;
811	}
812	return matrix;
813}
814
815/*---------------------------------------------------------------------*/
816
817function mapSets (sourceList,targetList)
818// source ID -> target ID (-1 means no correspondence)
819
820{
821	targetIndexing = {};
822	_d = Rows (targetList) * Columns (targetList);
823
824	for (_i = 0; _i < _d; _i += 1)
825	{
826		targetIndexing [targetList[_i]] = _i + 1;
827	}
828	_d = Rows (sourceList) * Columns (sourceList);
829	mapping 	  = {1,_d};
830	for (_i = 0; _i < _d; _i += 1)
831	{
832		mapping [_i] = targetIndexing[sourceList[_i]] - 1;
833	}
834
835	return mapping;
836}
837
838/*---------------------------------------------------------------------*/
839
840function mapStrings (sourceStr,targetStr)
841// source ID -> target ID (-1 means no correspondence)
842
843{
844	mapping 	  = {};
845	targetIndexing = {};
846	_d = Abs(targetStr);
847
848	for (_i = 0; _i < _d; _i += 1)
849	{
850		targetIndexing [targetStr[_i]] = _i + 1;
851	}
852	_d = Abs (sourceStr);
853	for (_i = 0; _i < _d; _i += 1)
854	{
855		mapping [_i] = targetIndexing[sourceStr[_i]] - 1;
856	}
857
858	return mapping;
859}
860/*---------------------------------------------------------------------*/
861
862function remapSequenceCoordinatesToReference (ref,seq)
863{
864	_seqLen	  = Abs(seq);
865	_coordMap = {_seqLen,1}["-1"];
866
867	_k				= (ref$"^\\-+");
868	_referenceSpan	= _k[1]+1;
869
870	for (_k = 0; _k < _referenceSpan; _k += 1) {
871		_coordMap[_k] = 0;
872	}
873
874	_qryCoord = _k;
875	_refCoord = 0;
876
877	while (_k < Abs(seq)) {
878		if (seq[_k] != "-") {
879			_coordMap[_qryCoord] = _refCoord;
880			_qryCoord += 1;
881		}
882		if (ref[_k] != "-") {
883			_refCoord += 1;
884		}
885		_k += 1;
886	}
887	return _coordMap;
888}
889
890/*---------------------------------------------------------------------*/
891
892function runModule (module,options,suppressStdout)
893{
894	if (suppressStdout)
895	{
896		_gfr = GLOBAL_FPRINTF_REDIRECT;
897		GLOBAL_FPRINTF_REDIRECT = "/dev/null";
898	}
899	LoadFunctionLibrary (module,options);
900	if (suppressStdout)
901	{
902		GLOBAL_FPRINTF_REDIRECT = _gfr;
903	}
904}
905
906/*---------------------------------------------------------------------*/
907
908function _formatTimeString (secondCount)
909{
910	_hours = secondCount $3600;
911	if (_hours < 10)
912	{
913		_timeString = "0"+_hours;
914	}
915	else
916	{
917		_timeString = ""+_hours;
918	}
919	_minutes = (secondCount%3600)$60;
920	if (_minutes < 10)
921	{
922		_timeString = _timeString+":0"+_minutes;
923	}
924	else
925	{
926		_timeString = _timeString+":"+_minutes;
927	}
928	_seconds = (secondCount%60);
929	if (_seconds < 10)
930	{
931		_timeString = _timeString+":0"+_seconds;
932	}
933	else
934	{
935		_timeString = _timeString+":"+_seconds;
936	}
937	return _timeString;
938}
939
940/*---------------------------------------------------------------------*/
941
942lfunction _constrainVariablesAndDescendants (variable) {
943    GetInformation (allVars, "^" + (variable&&6) + "\\..+$");
944    for (k = 0; k < Columns (allVars); k += 1) {
945        variableID    = allVars[k];
946        current_value = ^variableID;
947        ^variableID := current_value__;
948    }
949    return 0;
950}
951
952/*---------------------------------------------------------------------*/
953
954lfunction _unconstrainVariablesAndDescendants (variable) {
955    GetInformation (allVars, "^" + (variable&&6) + "\\..+$");
956    for (k = 0; k < Columns (allVars); k += 1) {
957        variableID    = allVars[k];
958        ClearConstraints (^variableID);
959    }
960    return 0;
961}
962