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