1LoadFunctionLibrary("../all-terms.bf"); 2LoadFunctionLibrary("parameters.bf"); 3LoadFunctionLibrary("../UtilityFunctions.bf"); 4LoadFunctionLibrary("model_functions.bf"); 5 6/** @module frequencies */ 7 8/** 9 * Sets model's equilibrium frequency estimator to Binary chatacter 2x1 estimator 10 * @name frequencies.empirical.binary 11 * @param {Dictionary} model 12 * @param {String} namespace - does nothing 13 * @param {DataSetFilter} datafilter - does nothing 14 * @returns {Dictionary} updated model 15 */ 16function frequencies.empirical.binary(model, namespace, datafilter) { 17 model = frequencies._aux.empirical.singlechar(model, namespace, datafilter); 18 model[terms.model.efv_estimate_name] = terms.frequencies.binary; 19 return model; 20} 21 22/** 23 * Sets model's equilibrium frequency estimator to equal frequencies 24 * @name frequencies.equal 25 * @param {Dictionary} model 26 * @param {String} namespace - does nothing 27 * @param {DataSetFilter} datafilter - does nothing 28 * @returns {Dictionary} updated model 29 */ 30function frequencies.equal(model, namespace, datafilter) { 31 32 33 __N = utility.Array1D(model[terms.alphabet]); 34 model[terms.efv_estimate] = { 35 __N, 36 1 37 }["1/__N"]; 38 model[terms.model.efv_estimate_name] = terms.frequencies.equal; 39 (model[terms.parameters])[terms.model.empirical] = 0; 40 return model; 41} 42 43/** 44 * Sets model's equilibrium frequency estimator to Nucleotide 4x1 estimator 45 * @name frequencies.empirical.nucleotide 46 * @param {Dictionary} model 47 * @param {String} namespace 48 * @param {DataSetFilter} datafilter 49 * @returns {Dictionary} updated model 50 */ 51function frequencies.empirical.nucleotide(model, namespace, datafilter) { 52 model = frequencies._aux.empirical.singlechar(model, namespace, datafilter); 53 model[terms.model.efv_estimate_name] = terms.frequencies._4x1; 54 (model[terms.parameters])[terms.model.empirical] = 3; 55 return model; 56} 57 58/** 59 * Compute equilibrium frequencies at run-time using Q inversion 60 * @name frequencies.empirical.nucleotide 61 * @param {Dictionary} model 62 * @param {String} namespace 63 * @param {DataSetFilter} datafilter 64 * @returns {Dictionary} updated model 65 */ 66 67lfunction frequencies.runtime.nucleotide(model, namespace, datafilter) { 68 69 /* 70 The run-time calculation depends on the current value of the rate matrix and 71 proceeds by 72 73 (1) setting the LAST column of the Q matrix to all 1s 74 (2) inverting the modified matrix 75 (3) reading off equilibrium frequencies from the last row of the inverse 76 77 This function needs to have the rate matrix instantiated before it can set 78 up the dependancies. 79 */ 80 81 82 return model; 83} 84 85/** 86 * This is an function that computes stationary frequencies of Markov process based on its 87 * rate matrix 88 * @param {Matrix} Q matrix 89*/ 90 91lfunction frequencies._aux.invert_model (Q) { 92 dim = Rows(Q); 93 for (i = 0; i < dim; i++) { 94 Q[i][dim-1] = 1; 95 } 96 return (Inverse (Q))[{{dim-1,0}}][{{dim-1,dim-1}}]; 97} 98 99 100/** 101 * Sets model's equilibrium frequency estimator to ML for binary data 102 * @name frequencies.empirical.binary 103 * @param {Dictionary} model 104 * @param {String} namespace 105 * @param {DataSetFilter} datafilter 106 * @returns {Dictionary} updated model 107 */ 108function frequencies.ML.binary (model, namespace, datafilter) { 109 model = frequencies._aux.empirical.singlechar(model, namespace, datafilter); 110 // define 2 frequency parameters 111 // initialize to empirical freqs 112 // add to model parameter manifest 113 frequencies.ML.binary.emp = model[terms.efv_estimate]; 114 model[terms.efv_estimate] = {2,1}; 115 frequencies.ML.binary.variables = {2,1}; 116 frequencies.ML.binary.scaler = namespace + ".frequency_scaler"; 117 118 utility.ForEachPair (model[terms.alphabet], "_index", "_letter", 119 ' 120 _idx = _index[1]; 121 _freq_parameter = namespace + ".equilibrium_frequency_of." + _letter; 122 frequencies.ML.binary.variables [_idx] = _freq_parameter; 123 (model[terms.efv_estimate]) [_idx] = _freq_parameter + "/" + frequencies.ML.binary.scaler; 124 parameters.DeclareGlobalWithRanges (_freq_parameter, frequencies.ML.binary.emp[_idx], 0, 1); 125 model.generic.AddGlobal (model, _freq_parameter, terms.characterFrequency (_letter)); 126 ' 127 ); 128 129 // constrain pi_0 + pi_1 = 1, 130 131 parameters.SetConstraint ( frequencies.ML.binary.scaler, Join ("+", frequencies.ML.binary.variables), "global"); 132 133 model[terms.model.efv_estimate_name] = terms.frequencies.ml; 134 (model[terms.parameters])[terms.model.empirical] = -1; // correct for the restricted sum 135 return model; 136} 137 138/** 139 * Sets model's equilibrium frequency estimator to protein 20x1 estimator 140 * @name frequencies.empirical.protein 141 * @param {Dictionary} model 142 * @param {String} namespace 143 * @param {DataSetFilter} datafilter 144 * @returns {Dictionary} updated model 145 */ 146function frequencies.empirical.protein (model, namespace, datafilter) { 147 model = frequencies._aux.empirical.singlechar(model, namespace, datafilter); 148 model[terms.model.efv_estimate_name] = terms.frequencies._20x1; 149 (model[terms.parameters])[terms.model.empirical] = 19; 150 return model; 151} 152 153/** 154 * Sets model's equilibrium frequency estimator to ML for protein data 155 * @name frequencies.empirical.protein 156 * @param {Dictionary} model 157 * @param {String} namespace 158 * @param {DataSetFilter} datafilter 159 * @returns {Dictionary} updated model 160 */ 161function frequencies.ML.protein (model, namespace, datafilter) { 162 return frequencies.mle (model, namespace, datafilter); 163} 164 165 166/** 167 * Multiplies nucleotide frequencies into the empirical codon estimate 168 * @name frequencies.empirical.protein 169 * @param {Dictionary} model 170 * @param {Matrix} __estimates 4x3 matrix of frequency estimates 171 * @returns {Dictionary} updated model 172 */ 173function frequencies.codon.multiply_in_frequencies (model, __estimates) { 174 175 __dimension = model.Dimension(model); 176 __alphabet = model[^"terms.alphabet"]; 177 178 if (Type (model[terms.model.rate_matrix]) == "AssociativeList") { // handle mixture models 179 __components = Rows (model[terms.model.rate_matrix]); 180 __component_count = utility.Array1D (__components); 181 182 for (_rowChar = 0; _rowChar < __dimension; _rowChar += 1) { 183 for (_colChar = _rowChar + 1; _colChar < __dimension; _colChar += 1) { 184 185 //__diff = models.codon.diff(__alphabet[_rowChar], __alphabet[_colChar]); 186 for (__component_id = 0; __component_id < __component_count; __component_id += 1) { 187 if (Abs (((model[terms.model.rate_matrix])[__components[__component_id]]) [_colChar][_rowChar])) { 188 __diff = models.codon.diff.complete (__alphabet[_rowChar], __alphabet[_colChar]); 189 for (__diff_id = 0; __diff_id < Abs (__diff); __diff_id += 1) { 190 ((model[terms.model.rate_matrix])[__components[__component_id]]) [_rowChar][_colChar] += "*" + (__estimates[(__diff[__diff_id])["to"]])[(__diff[__diff_id])["position"]]; 191 ((model[terms.model.rate_matrix])[__components[__component_id]]) [_colChar][_rowChar] += "*" + (__estimates[(__diff[__diff_id])["from"]])[(__diff[__diff_id])["position"]]; 192 } 193 } 194 } 195 } 196 } 197 } else { 198 for (_rowChar = 0; _rowChar < __dimension; _rowChar += 1) { 199 for (_colChar = _rowChar + 1; _colChar < __dimension; _colChar += 1) { 200 201 if (Abs ((model[terms.model.rate_matrix])[_colChar][_rowChar])) { 202 __diff = models.codon.diff.complete (__alphabet[_rowChar], __alphabet[_colChar]); 203 for (__diff_id, __diff_value; in; __diff) { 204 (model[terms.model.rate_matrix])[_rowChar][_colChar] += "*" + (__estimates[__diff_value["to"]])[__diff_value["position"]]; 205 (model[terms.model.rate_matrix])[_colChar][_rowChar] += "*" + (__estimates[__diff_value["from"]])[__diff_value["position"]]; 206 } 207 } 208 } 209 } 210 } 211 return model; 212} 213 214 215 216lfunction frequencies._aux.validate (model) { 217 218 __dimension = model.Dimension(model); 219 __alphabet = model[^"terms.alphabet"]; 220 221 if (Type (model[^"terms.model.rate_matrix"]) == "AssociativeList") { 222 utility.ForEach (model[^"terms.model.rate_matrix"], "_q_matrix_", ' 223 assert(Type(_q_matrix_) == "Matrix" && Rows(_q_matrix_) == __dimension && Columns(_q_matrix_) == __dimension, 224 "`^'terms.model.rate_matrix'` must be defined prior to calling frequency estimators") 225 '); 226 } else { 227 assert(Type(model[^"terms.model.rate_matrix"]) == "Matrix" && Rows(model[^"terms.model.rate_matrix"]) == __dimension && Columns(model[^"terms.model.rate_matrix"]) == __dimension, 228 "`^'terms.model.rate_matrix'` must be defined prior to calling frequency estimators"); 229 } 230 231} 232 233lfunction frequencies.empirical.codon_from_nuc (model, nuc_dict) { 234 235 236 result = {model.Dimension (model),1}; 237 corrector = 1; 238 utility.ForEach (model[^"terms.stop_codons"], "codon", 239 ' 240 `&corrector` += - (`&nuc_dict`[codon[0]])[0] * (`&nuc_dict`[codon[1]])[1] * (`&nuc_dict`[codon[2]])[2]; 241 '); 242 243 utility.ForEachPair (model[^"terms.alphabet"], "index", "codon", 244 ' 245 (`&result`)[index[1]] = (`&nuc_dict`[codon[0]])[0] * (`&nuc_dict`[codon[1]])[1] * (`&nuc_dict`[codon[2]])[2] / `&corrector`; 246 '); 247 248 return (result); 249} 250/** 251 * @name frequencies.empirical.F3x4 252 * @param {Dictionary} model 253 * @param {String} namespace 254 * @param {DataSetFilter} datafilter 255 * @returns {Dictionary} updated model 256 */ 257 258lfunction frequencies.empirical.F3x4(model, namespace, datafilter) { 259 260 frequencies._aux.validate (model); 261 262 if (None != datafilter) { 263 utility.ToggleEnvVariable("COUNT_GAPS_IN_FREQUENCIES", FALSE); 264 __f = frequencies._aux.empirical.collect_data(datafilter, 3, 1, 1); 265 utility.ToggleEnvVariable("COUNT_GAPS_IN_FREQUENCIES", None); 266 } else { 267 __f = model [^"terms.efv_estimate"]; 268 } 269 270 271 __alphabet = model[^"terms.bases"]; 272 nuc_dict = {}; 273 utility.ForEachPair (__alphabet, "index","base", 274 ' 275 `&nuc_dict`[base] = `&__f`[index[1]][-1]; 276 '); 277 278 frequencies.codon.multiply_in_frequencies (model, nuc_dict); 279 model[^"terms.model.efv_estimate_name"] = ^"terms.frequencies.F3x4"; 280 model[^"terms.efv_estimate"] = frequencies.empirical.codon_from_nuc (model, nuc_dict); 281 (model[^"terms.parameters"])[^"terms.model.empirical"] = 9; 282 283 return model; 284} 285 286/** 287 * @name frequencies.empirical.F1x4 288 * @param {Dictionary} model 289 * @param {String} namespace 290 * @param {DataSetFilter} datafilter 291 * @returns {Dictionary} updated model 292 */ 293 294lfunction frequencies.empirical.F1x4(model, namespace, datafilter) { 295 296 frequencies._aux.validate (model); 297 298 if (None != datafilter) { 299 utility.ToggleEnvVariable("COUNT_GAPS_IN_FREQUENCIES", FALSE); 300 __f1 = frequencies._aux.empirical.collect_data(datafilter, 1, 1, 1); 301 __f = {4,3} ["__f1[_MATRIX_ELEMENT_ROW_][0]"]; 302 utility.ToggleEnvVariable("COUNT_GAPS_IN_FREQUENCIES", None); 303 } else { 304 __f = model [^"terms.efv_estimate"]; 305 } 306 307 __alphabet = model[^"terms.bases"]; 308 nuc_dict = {}; 309 utility.ForEachPair (__alphabet, "index","base", 310 ' 311 `&nuc_dict`[base] = `&__f`[index[1]][-1]; 312 '); 313 314 frequencies.codon.multiply_in_frequencies (model, nuc_dict); 315 model[^"terms.model.efv_estimate_name"] = ^"terms.frequencies.F1x4"; 316 model[^"terms.efv_estimate"] = frequencies.empirical.codon_from_nuc (model, nuc_dict); 317 (model[^"terms.parameters"])[^"terms.model.empirical"] = 3; 318 return model; 319} 320 321 322/** 323 * @name frequencies.empirical.corrected.CF3x4 324 * @param {Dictionary} model 325 * @param {String} namespace 326 * @param {DataSetFilter} datafilter 327 * @returns {Dictionary} updated model 328 */ 329lfunction frequencies.empirical.corrected.CF3x4(model, namespace, datafilter) { 330 331 frequencies._aux.validate (model); 332 333 if (None != datafilter) { 334 utility.ToggleEnvVariable("COUNT_GAPS_IN_FREQUENCIES", 0); 335 __f = frequencies._aux.empirical.collect_data(datafilter, 3, 1, 1); 336 utility.ToggleEnvVariable("COUNT_GAPS_IN_FREQUENCIES", None); 337 } else { 338 __f = model [^"terms.efv_estimate"]; 339 } 340 341 //TODO 342 __alphabet = model[^"terms.alphabet"]; 343 __estimates = frequencies._aux.CF3x4(__f, model[^"terms.bases"], __alphabet, model[^"terms.stop_codons"]); 344 345 model[^"terms.efv_estimate"] = __estimates[^"terms.codons"]; 346 347 frequencies.codon.multiply_in_frequencies (model, __estimates[^"terms.bases"]); 348 349 model[^"terms.model.efv_estimate_name"] = ^"terms.frequencies.CF3x4"; 350 (model[^"terms.parameters"])[^"terms.model.empirical"] = 9; 351 return model; 352} 353 354/** 355 * To be implemented 356 * @name frequencies.mle 357 * @param model 358 * @param namespace 359 * @param datafilter 360 */ 361lfunction frequencies.mle(model, namespace, datafilter) { 362 frequencies._aux.empirical.singlechar(model, namespace, datafilter); // gather empirical frequencies 363 364 alphabet = model[utility.getGlobalValue ("terms.alphabet")]; 365 alphabet_size = utility.Array1D (alphabet); 366 empirical = model[utility.getGlobalValue ("terms.efv_estimate")]; 367 368 model[utility.getGlobalValue ("terms.efv_estimate")] = {alphabet_size,1}; 369 370 ML.variables = {alphabet_size,1}; 371 ML.scaler = namespace + ".frequency_scaler"; 372 373 utility.ForEachPair (model[utility.getGlobalValue ("terms.alphabet")], "_index", "_letter", 374 ' 375 _idx = _index[1]; 376 _freq_parameter = `&namespace` + ".equilibrium_frequency_of." + _letter; 377 (`&ML.variables`) [_idx] = _freq_parameter; 378 (`&model`[terms.efv_estimate]) [_idx] = _freq_parameter + "/" + `&ML.scaler`; 379 parameters.DeclareGlobalWithRanges (_freq_parameter, `&empirical`[_idx], 0, 1); 380 model.generic.AddGlobal (`&model`, _freq_parameter, terms.characterFrequency (_letter)); 381 ' 382 ); 383 384 // constrain pi_A + ... + pi_A = 1, 385 386 parameters.SetConstraint ( ML.scaler, Join ("+", ML.variables), "global"); 387 388 model[utility.getGlobalValue("terms.model.efv_estimate_name")] = utility.getGlobalValue ("terms.frequencies.MLE"); 389 (model[utility.getGlobalValue("terms.parameters")])[utility.getGlobalValue("terms.model.empirical")] = -1; // correct for the restricted sum 390 391 return model; 392 393} 394 395//--- AUX FUNCTIONS FROM THIS POINT ON ---// 396 397/** 398 * Returns entire character count (# of sites * # of species) 399 * @name frequencies._aux.empirical.character_count 400 * @private 401 * @param {DataSetFilter} datafilter - the datafilter containing site and species information 402 */ 403function frequencies._aux.empirical.character_count(datafilter) { 404 return Eval("`datafilter`.sites * `datafilter`.species"); 405} 406 407/** 408 * Collects frequency data from dataset 409 * @name frequencies._aux.empirical.collect_data 410 * @private 411 * @param {DataSetFilter} datafilter 412 * @param {Number} unit 413 * @param {Number} stride 414 * @param {Number} position_specific 415 * @returns {Matrix} frequency data 416 */ 417lfunction frequencies._aux.empirical.collect_data(datafilter, unit, stride, position_specific) { 418 419 assert(Type(datafilter) == "Matrix" || Type(datafilter) == "AssociativeList" || Type(datafilter) == "String", 420 "datafilter must be a String/Matrix/Associative List in call to frequencies._aux.empirical.collect_data, instead had a " + Type(datafilter) 421 ); 422 if (Type(datafilter) == "String") { 423 HarvestFrequencies(__f, ^ datafilter, unit, stride, position_specific); 424 } else { 425 site_count = 0; 426 dim = utility.Array1D(datafilter); 427 428 for (i = 0; i < dim; i += 1) { 429 430 HarvestFrequencies(__f, ^ (datafilter[i]), unit, stride, position_specific); 431 local_sites = frequencies._aux.empirical.character_count(datafilter[i]); 432 433 if (i) { 434 __f_composite += __f * local_sites; 435 } else { 436 __f_composite = __f * local_sites; 437 438 } 439 site_count += local_sites; 440 } 441 return __f_composite * (1 / site_count); 442 } 443 444 445 return __f; 446} 447 448/** 449 * Updates model's equilibrium frequency estimator to not count gaps in frequencies 450 * @name frequencies._aux.empirical.singlechar 451 * @private 452 * @param {Dictionary} model - the model to update the 453 * @param {String} namespace - does nothing 454 * @param {DataSetFilter} datafilter - the datasetfilter to collect data from 455 * @returns {Dictionary} updated model 456 */ 457function frequencies._aux.empirical.singlechar(model, namespace, datafilter) { 458 utility.ToggleEnvVariable("COUNT_GAPS_IN_FREQUENCIES", 0); 459 __f = frequencies._aux.empirical.collect_data(datafilter, 1, 1, 1); 460 utility.ToggleEnvVariable("COUNT_GAPS_IN_FREQUENCIES", None); 461 model[terms.efv_estimate] = __f; 462 return model; 463} 464 465/** 466 * @name frequencies._aux.CF3x4 467 * @private 468 * @param observed_3x4 469 * @param base_alphabet 470 * @param sense_codons 471 * @param stop_codons 472 * @returns {Dictionary} codons and bases 473 */ 474function frequencies._aux.CF3x4(observed_3x4, base_alphabet, sense_codons, stop_codons) { 475 476 frequencies._aux.CF3x4.p = {}; 477 478 frequencies._aux.CF3x4.args = {}; 479 480 for (frequencies._aux.CF3x4.k = 0; frequencies._aux.CF3x4.k < 3; frequencies._aux.CF3x4.k += 1) { 481 frequencies._aux.CF3x4.p[frequencies._aux.CF3x4.k] = parameters.GenerateSequentialNames("frequencies._aux.CF3x4.p" + frequencies._aux.CF3x4.k, 3, None); 482 parameters.DeclareGlobal(frequencies._aux.CF3x4.p[frequencies._aux.CF3x4.k], None); 483 parameters.SetRange(frequencies._aux.CF3x4.p[frequencies._aux.CF3x4.k], terms.range01); 484 frequencies._aux.CF3x4.args + (Join(",", frequencies._aux.CF3x4.p[frequencies._aux.CF3x4.k])); 485 } 486 487 frequencies._aux.CF3x4.args = Join(",", frequencies._aux.CF3x4.args); 488 489 frequencies._aux.CF3x4.n = {}; 490 491 for (frequencies._aux.CF3x4.k = 0; frequencies._aux.CF3x4.k < 3; frequencies._aux.CF3x4.k += 1) { 492 frequencies._aux.CF3x4.n[frequencies._aux.CF3x4.k] = parameters.GenerateAttributedNames("frequencies._aux.CF3x4.n" + frequencies._aux.CF3x4.k, base_alphabet, None); 493 parameters.SetConstraint(frequencies._aux.CF3x4.n[frequencies._aux.CF3x4.k], 494 parameters.helper.stick_breaking(frequencies._aux.CF3x4.p[frequencies._aux.CF3x4.k], observed_3x4[-1][frequencies._aux.CF3x4.k]), 495 utility.getGlobalValue("terms.global")); 496 } 497 498 499 frequencies._aux.stop_count = Columns(stop_codons); 500 501 frequencies._aux.CF3x4.stop_correction = {}; 502 503 for (frequencies._aux.CF3x4.i = 0; frequencies._aux.CF3x4.i < Columns(base_alphabet); frequencies._aux.CF3x4.i += 1) { 504 frequencies._aux.CF3x4.stop_correction[base_alphabet[frequencies._aux.CF3x4.i]] = { 505 { 506 "", 507 "", 508 "" 509 } 510 }; 511 } 512 513 frequencies._aux.CF3x4.denominator = "1"; 514 515 for (frequencies._aux.CF3x4.i = 0; frequencies._aux.CF3x4.i < frequencies._aux.stop_count; frequencies._aux.CF3x4.i += 1) { 516 frequencies._aux.CF3x4.sc = stop_codons[frequencies._aux.CF3x4.i]; 517 518 (frequencies._aux.CF3x4.stop_correction[frequencies._aux.CF3x4.sc[0]])[0] += 519 "-frequencies._aux.CF3x4.n1_" + frequencies._aux.CF3x4.sc[1] + 520 "*frequencies._aux.CF3x4.n2_" + frequencies._aux.CF3x4.sc[2]; 521 522 (frequencies._aux.CF3x4.stop_correction[frequencies._aux.CF3x4.sc[1]])[1] += 523 "-frequencies._aux.CF3x4.n0_" + frequencies._aux.CF3x4.sc[0] + 524 "*frequencies._aux.CF3x4.n2_" + frequencies._aux.CF3x4.sc[2]; 525 526 (frequencies._aux.CF3x4.stop_correction[frequencies._aux.CF3x4.sc[2]])[2] += 527 "-frequencies._aux.CF3x4.n0_" + frequencies._aux.CF3x4.sc[0] + 528 "*frequencies._aux.CF3x4.n1_" + frequencies._aux.CF3x4.sc[1]; 529 530 frequencies._aux.CF3x4.denominator += "-frequencies._aux.CF3x4.n0_" + frequencies._aux.CF3x4.sc[0] + 531 "*frequencies._aux.CF3x4.n1_" + frequencies._aux.CF3x4.sc[1] + 532 "*frequencies._aux.CF3x4.n2_" + frequencies._aux.CF3x4.sc[2]; 533 } 534 535 536 parameters.SetConstraint("frequencies._aux.CF3x4.denominator", frequencies._aux.CF3x4.denominator, utility.getGlobalValue("terms.global")); 537 538 frequencies._aux.N = { 539 Columns(base_alphabet), 540 3 541 }; 542 frequencies._aux.res = {}; 543 frequencies._aux.codons = { 544 Columns(sense_codons), 545 1 546 }; 547 548 for (frequencies._aux.CF3x4.i = 0; frequencies._aux.CF3x4.i < Columns(sense_codons); frequencies._aux.CF3x4.i += 1) { 549 frequencies._aux.CF3x4.sc = { 550 3, 551 1 552 }; 553 for (frequencies._aux.CF3x4.pos = 0; frequencies._aux.CF3x4.pos < 3; frequencies._aux.CF3x4.pos += 1) { 554 frequencies._aux.CF3x4.sc[frequencies._aux.CF3x4.pos] = "frequencies._aux.CF3x4.n" + frequencies._aux.CF3x4.pos + "_" + (sense_codons[frequencies._aux.CF3x4.i])[frequencies._aux.CF3x4.pos]; 555 } 556 ExecuteCommands("frequencies._aux.codons[frequencies._aux.CF3x4.i] := " + Join("*", frequencies._aux.CF3x4.sc) + "/frequencies._aux.CF3x4.denominator"); 557 } 558 559 for (frequencies._aux.CF3x4.i = 0; frequencies._aux.CF3x4.i < Columns(base_alphabet); frequencies._aux.CF3x4.i += 1) { 560 frequencies._aux.CF3x4.n = base_alphabet[frequencies._aux.CF3x4.i]; 561 frequencies._aux.res[frequencies._aux.CF3x4.n] = { 562 3, 563 1 564 }; 565 for (frequencies._aux.CF3x4.pos = 0; frequencies._aux.CF3x4.pos < 3; frequencies._aux.CF3x4.pos += 1) { 566 567 frequencies._aux.CF3x4.sc = (frequencies._aux.CF3x4.stop_correction[frequencies._aux.CF3x4.n])[frequencies._aux.CF3x4.pos]; 568 569 if (Abs(frequencies._aux.CF3x4.sc)) { 570 frequencies._aux.CF3x4.sc = "*(1" + frequencies._aux.CF3x4.sc + ")"; 571 } 572 573 ExecuteCommands("frequencies._aux.N[" + frequencies._aux.CF3x4.i + "][" + frequencies._aux.CF3x4.pos + "] := frequencies._aux.CF3x4.n" + frequencies._aux.CF3x4.pos + "_" + frequencies._aux.CF3x4.n + frequencies._aux.CF3x4.sc + "/frequencies._aux.CF3x4.denominator"); 574 575 ExecuteCommands("(frequencies._aux.res[frequencies._aux.CF3x4.n])[frequencies._aux.CF3x4.pos] := frequencies._aux.CF3x4.n" + frequencies._aux.CF3x4.pos + "_" + frequencies._aux.CF3x4.n); 576 } 577 } 578 579 580 581 ExecuteCommands("Optimize (frequencies._aux.CF3x4.p, frequencies._aux._CF3x4_minimizer( " + 582 frequencies._aux.CF3x4.args + "))"); 583 584 return { 585 utility.getGlobalValue("terms.codons"): Eval("frequencies._aux.codons"), 586 utility.getGlobalValue("terms.bases"): utility.Map (frequencies._aux.res,"_value_", "Eval(_value_)") // this is an in-place shallow copy 587 }; 588} 589 590/** 591 * @name frequencies._aux._CF3x4_minimizer 592 * @private 593 * @param p11 594 * @param p12 595 * @param p13 596 * @param p21 597 * @param p22 598 * @param p23 599 * @param p31 600 * @param p32 601 * @param p33 602 */ 603function frequencies._aux._CF3x4_minimizer(p11, p12, p13, p21, p22, p23, p31, p32, p33) { 604 frequencies._aux._CF3x4_minimizer.error = frequencies._aux.N - observed_3x4; 605 return -(+frequencies._aux._CF3x4_minimizer.error$frequencies._aux._CF3x4_minimizer.error); 606} 607