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