1skipCodeSelectionStep 		= 0;
2
3
4LoadFunctionLibrary("chooseGeneticCode");
5LoadFunctionLibrary("GrabBag");
6LoadFunctionLibrary("dSdNTreeTools");
7LoadFunctionLibrary("CF3x4");
8LoadFunctionLibrary("BranchSiteTemplate");
9LoadFunctionLibrary("DescriptiveStatistics");
10LoadFunctionLibrary("TreeTools");
11
12
13DataSet 			ds 				= ReadDataFile(PROMPT_FOR_FILE);
14DataSetFilter 		dsf 			= CreateFilter(ds,3,"","",GeneticCodeExclusions);
15
16HarvestFrequencies	(nuc3, dsf, 3, 1, 1);
17nucCF			  = CF3x4	(nuc3, GeneticCodeExclusions);
18
19rate_class_count = prompt_for_a_value ("How many rate classes", 3, 2, 8, 1);
20fprintf (stdout, "\nUsing ", rate_class_count, " rate classes\n");
21
22PopulateModelMatrix			  ("MGMatrixLocal",  nucCF, "syn", "", "nonsyn");
23codon3x4					= BuildCodonFrequencies (nucCF);
24Model		MGL				= (MGMatrixLocal, codon3x4, 0);
25
26_DO_TREE_REBALANCE_ = 0;
27LoadFunctionLibrary			  ("queryTree");
28
29SetDialogPrompt ("Save analysis results to");
30fprintf (PROMPT_FOR_FILE, CLEAR_FILE, KEEP_OPEN,"Model\tLogL\tNP\tBIC\tTree String");
31csvFilePath = LAST_FILE_PATH;
32
33
34// ---- PHASE 0 : LOCAL (no omega variation model) FIT. -----
35
36fprintf 					  (stdout, "[BS-REL PHASE 0] Fitting the local MG94 (no site-to-site variation) to obtain initial parameter estimates\n");
37
38
39LikelihoodFunction	base_LF	 = (dsf, givenTree);
40
41/*
42ExecuteAFile (csvFilePath + ".mglocal.fit");
43res_base = {2,2};
44*/
45
46Optimize					  (res_base,base_LF);
47writeTheLF (".mglocal.fit", "base_LF");
48baseParameters                = 9;
49localLL						 = res_base[1][0];
50localParams					 = res_base[1][1] + baseParameters;
51
52totalBranchCount			 = BranchCount(givenTree) + TipCount (givenTree);
53
54pValueByBranch				  = {totalBranchCount,10};
55bNames						  = BranchName (givenTree, -1);
56
57for (k = 0; k < totalBranchCount; k = k+1) {
58	srate  = Eval ("givenTree." + bNames[k] + ".syn");
59	nsrate = Eval ("givenTree." + bNames[k] + ".nonsyn");
60	if (srate > 0) {
61		pValueByBranch [k][0] = Min (10, nsrate/srate);
62	} else {
63		pValueByBranch [k][0] = 10;
64	}
65}
66
67sample_size                    = dsf.sites * dsf.species;
68omegaStats					 = GatherDescriptiveStats (pValueByBranch[-1][0]);
69fprintf						 (stdout, "Local omega model (no rate variation) ");
70
71
72localBIC                    = BIC (localLL, localParams, sample_size);
73PrintDescriptiveStats		 ("Branch omega values", omegaStats);
74
75for (rcc = 1; rcc <= rate_class_count; rcc += 1) {
76    PopulateModelMatrix			  ("MGMatrix"+rcc,  nucCF, "t", "omega"+rcc, "");
77    ExecuteCommands               ("global omegaG" + rcc + " = " + (0.3*rcc - 0.2));
78    PopulateModelMatrix			  ("MGMatrixG"+rcc,  nucCF, "t", "omegaG"+rcc, "");
79}
80
81// ---- PHASE 1 : FULL GLOBAL MODEL FIT ---- //
82
83generate_bs_rel_gdd_freqs (rate_class_count, "global_freqsG", "model_mixingG", "PauxG", "MGMatrixG",1);
84generate_bs_rel_gdd_freqs (rate_class_count, "global_freqs",  "model_mixing",  "Paux", "MGMatrix",0);
85
86sortedLocalDNDS              = (pValueByBranch[-1][0])%0;
87
88for (rcc = 1; rcc <= rate_class_count; rcc += 1) {
89    ExecuteCommands ("omegaG" + rcc + "=" +  sortedLocalDNDS[Max(0,totalBranchCount*rcc$rate_class_count-1)]);
90}
91
92
93ExecuteCommands ("Model 		MGG	=(\"`model_mixingG`\",codon3x4,EXPLICIT_FORM_MATRIX_EXPONENTIAL);");
94Tree						   mixtureTreeG = treeString;
95
96ExecuteCommands ("Model 		MGLocal=(\"`model_mixing`\",codon3x4,EXPLICIT_FORM_MATRIX_EXPONENTIAL);");
97Tree						   mixtureTree = treeString;
98
99ReplicateConstraint 		  ("this1.?.t:=this2.?.syn",mixtureTreeG,givenTree);
100ReplicateConstraint 		  ("this1.?.t:=this2.?.syn",mixtureTree,givenTree);
101ClearConstraints			  (mixtureTree);
102ClearConstraints			  (mixtureTreeG);
103
104
105ASSUME_REVERSIBLE_MODELS	  = 1;
106OPTIMIZATION_METHOD           = 4;
107
108BIC_scores                    = {};
109BranchLengthEstimates         = {};
110
111mixTreeAVL                    = mixtureTreeG ^ 0;
112
113LikelihoodFunction three_LFG   = (dsf,mixtureTreeG);
114fprintf 					  (stdout, "[BS-REL PHASE 1] Fitting a GLOBAL branch-site matrix mixture\n");
115
116
117blG = getBranchLengthExpression (1);
118removeT = "/(" + (blG ^{{"\\*t",""}}) + ")";
119
120fprintf (stdout, base_LF);
121
122blLocal = {};
123for (k = 0; k < totalBranchCount; k += 1) {
124    blLocal [bNames[k]] = Eval("BranchLength(givenTree,\""+bNames[k]+"\")")*3;
125    ExecuteCommands ("mixtureTreeG." + bNames[k] + ".t := " + blLocal [bNames[k]] + removeT );
126}
127
128
129USE_LAST_RESULTS			= 1;
130SHORT_MPI_RETURN            = 1;
131VERBOSITY_LEVEL				= 0;
132
133OPTIMIZATION_PRECISION      = 0.1;
134
135Optimize					  (res_three_LF_global,three_LFG);
136
137fprintf (stdout, "\n\nApproximation phase 1 (use MG-local branch lengths): ", res_three_LF_global[1][0], "\n");
138fprintf     (stdout,"Tree branch lengths:\n",PostOrderAVL2StringDistances(mixTreeAVL,getBranchLengths ("mixtureTreeG", 1)));
139
140OPTIMIZATION_PRECISION      = 0.01;
141
142for (k = 0; k < totalBranchCount; k += 1) {
143    ExecuteCommands ("mixtureTreeG." + bNames[k] + ".t = mixtureTreeG." + bNames[k] + ".t");
144    bl = blLocal [bNames[k]] ;
145    ExecuteCommands ("FindRoot (z,`blG`-"+bl+",t,0,10000);");
146    ExecuteCommands ("mixtureTreeG." + bNames[k] + ".t :< " + 10*z);
147
148}
149
150Optimize					  (res_three_LF_global,three_LFG);
151fprintf (stdout, "\n\nApproximation phase 2 (maximum branch lengths are limited to 10x those of the null model): ", res_three_LF_global[1][0], "\n");
152fprintf     (stdout,"Tree branch lengths:\n",PostOrderAVL2StringDistances(mixTreeAVL,getBranchLengths ("mixtureTreeG", 1)));
153
154OPTIMIZATION_PRECISION      = 0.001;
155
156for (k = 0; k < totalBranchCount; k += 1) {
157    ExecuteCommands ("mixtureTreeG." + bNames[k] + ".t :< 10000");
158}
159
160Optimize					  (res_three_LF_global,three_LFG);
161
162//ExecuteAFile (csvFilePath + ".relglobal.fit");
163
164
165blStashByName = {};
166for (k = 0; k < totalBranchCount; k += 1) {
167    bl = blLocal [bNames[k]] ;
168    ExecuteCommands ("FindRoot (z,`blG`-"+bl+",t,0,10000);");
169    blStashByName[bNames[k]] = z;
170    //fprintf (stdout, "\n", bNames[k], " : ", z, " (", bl/3, ")\n");
171}
172
173writeTheLF (".relglobal.fit", "three_LFG");
174
175
176fprintf (stdout, "Global model fit:");
177
178BIC_scores ["global model"] = BIC (res_three_LF_global[1][0], res_three_LF_global[1][1]+baseParameters, sample_size);
179BranchLengthEstimates ["global model"] = getBranchLengths ("mixtureTreeG", 1);
180
181fprintf     (csvFilePath, "\nGlobal\t", res_three_LF_global[1][0],"\t", res_three_LF_global[1][1]+baseParameters, "\t", BIC_scores ["global model"], "\t", PostOrderAVL2StringDistances(mixTreeAVL,BranchLengthEstimates ["global model"]));
182
183fprintf (stdout, "\nInferred this global omega distribution: ");
184reportOmegaDistro ("");
185fprintf     (stdout,"Tree branch lengths:\n",PostOrderAVL2StringDistances(mixTreeAVL,BranchLengthEstimates ["global model"]));
186globalState = saveLF ("three_LFG");
187
188
189for (k = 0; k < totalBranchCount; k += 1) {
190    constrainABranch (bNames[k], 0);
191}
192
193// ---- PHASE 2 : LOCAL MODEL FITS ---- //
194
195fprintf (stdout, "\n[BS-REL PHASE 2] Fitting a GLOBAL branch-site matrix mixture with a SINGLE unconstrained branch\n");
196
197OPTIMIZATION_METHOD = 0;
198
199LikelihoodFunction three_LF   = (dsf,mixtureTree);
200
201branchValues = {};
202
203if (MPI_NODE_COUNT > 1) {
204    MPI_NODE_STATE = {MPI_NODE_COUNT-1,1};
205    MPI_NODE_STATE[0] = "";
206}
207
208
209for (k = 0; k < totalBranchCount; k+=1) {
210    thisBranchName = bNames[k];
211    fprintf (stdout, "\n[BS-REL PHASE 2. Branch '", thisBranchName, "']\n");
212    globalState["restoreLF"][""];
213    if (k) {
214        if (MPI_NODE_COUNT > 1) {
215            for (nodeID = 0; nodeID < totalBranchCount; nodeID += 1) {
216                constrainABranch (bNames[nodeID],0);
217            }
218        } else {
219            constrainABranch (bNames[k-1],0);
220        }
221    }
222    unConstrainABranch (thisBranchName);
223
224    if (MPI_NODE_COUNT > 1) {
225        for (nodeID = 0; nodeID < MPI_NODE_COUNT-1; nodeID += 1) {
226            if (Abs(MPI_NODE_STATE[nodeID]) == 0) {
227                MPISend (nodeID+1, three_LF);
228                fprintf (stdout, "\n[SENT TO NODE ", nodeID, "]\n");
229                MPI_NODE_STATE [nodeID] = thisBranchName;
230                break;
231            }
232        }
233        if (nodeID == MPI_NODE_COUNT-1) {
234            processABranch (thisBranchName, 1);
235        }
236    } else {
237        OPTIMIZATION_PRECISION = 0.05;
238        maxBL = blStashByName[thisBranchName]*10;
239        t = maxBL;
240        fprintf (stdout, thisBranchName, ": maximum t value: ", maxBL, ". Maximum achieved branch length: ", Eval (blG), "\n");
241        constrainABranch (thisBranchName,maxBL);
242        Optimize (localBranchRes, three_LF);
243        fprintf (stdout, "Approximation phase 1 (constrained branch length): ", localBranchRes[1][0], "\n");
244        constrainABranch (thisBranchName,-1);
245        OPTIMIZATION_PRECISION = 0.001;
246        Optimize (localBranchRes, three_LF);
247        processABranch (thisBranchName,0);
248    }
249}
250
251
252leftOver = 0;
253for (nodeID = 0; nodeID < MPI_NODE_COUNT-1; nodeID += 1) {
254    leftOver += Abs(MPI_NODE_STATE[nodeID])>0;
255}
256
257//fprintf (stdout, MPI_NODE_STATE, "\n", leftOver, "\n");
258
259
260for (nodeID = 0; nodeID < leftOver; nodeID += 1) {
261    processABranch ("", 0);
262}
263
264branchValues ["restoreLF"][""];
265
266fprintf (stdout, "\n[BS-REL PHASE 3] Fitting a LOCAL branch-site matrix mixture model\n");
267Optimize (res_local, three_LF);
268
269writeTheLF (".local.fit", "three_LF");
270
271fprintf (stdout, "Local model fit:");
272BIC_scores ["local model"] = BIC (res_local[1][0], res_local[1][1]+baseParameters, sample_size);
273BranchLengthEstimates ["local model"] = getBranchLengths ("mixtureTree", 0);
274fprintf     (stdout, "\nTree: ", PostOrderAVL2StringDistances(mixTreeAVL,BranchLengthEstimates ["local model"]), "\n");
275fprintf     (csvFilePath, "\nLocal\t", res_local[1][0],"\t", res_local[1][1]+baseParameters, "\t", BIC_scores ["local model"], "\t", PostOrderAVL2StringDistances(mixTreeAVL,BranchLengthEstimates ["local model"]), CLOSE_FILE);
276
277min_BIC = 1e100;
278
279function compute_min_BIC (key, value) {
280    min_BIC = Min (min_BIC, value);
281    return 0;
282}
283
284function akaike_weights_pass1 (key, value) {
285    akaike_scores [key] = Exp (0.5*(min_BIC-value));
286    return 0;
287}
288
289
290function akaike_weights_pass2 (key, value) {
291    BIC_scores [key] = value/norm_factor;
292    return 0;
293}
294
295fprintf (stdout, "\n[BS-REL PHASE 4] Generating a model averaged branch lengths\n");
296akaike_scores = {};
297BIC_scores ["compute_min_BIC"][""];
298BIC_scores ["akaike_weights_pass1"][""];
299norm_factor = +akaike_scores;
300akaike_scores ["akaike_weights_pass2"][""];
301
302reweighted_branch_lengths = {};
303
304model_names = Rows(BranchLengthEstimates);
305for (k = 0; k < Abs (BranchLengthEstimates); k+=1) {
306    key = model_names [k];
307    bl_avl = BranchLengthEstimates[key];
308    bnames = Rows (bl_avl);
309
310    for (b = 0; b < Abs (bl_avl); b += 1) {
311        reweighted_branch_lengths[bnames[b]] += BIC_scores[key]*bl_avl[bnames[b]];
312    }
313}
314
315fprintf (stdout,  PostOrderAVL2StringDistances(mixTreeAVL,reweighted_branch_lengths), "\n");
316
317//------------------------------------------------------------------------------------------------------------------------
318function processABranch (thisBranchName, doSend) {
319    if (MPI_NODE_COUNT > 1) {
320        MPIReceive (-1,fromNode,resStr);
321	    prevBranch = MPI_NODE_STATE[fromNode-1];
322        fprintf (stdout, "\n[RECEIVED " + prevBranch +" FROM NODE ", (fromNode-1), "]\n");
323        if (doSend) {
324            MPISend (fromNode, three_LF);
325            MPI_NODE_STATE [fromNode-1] = thisBranchName;
326            fprintf (stdout, "\n[SENT " + thisBranchName +" TO NODE ", (fromNode-1), "]\n");
327        } else {
328            MPI_NODE_STATE [fromNode-1] = "";
329        }
330        thisBranchName = prevBranch;
331        ExecuteCommands (resStr);
332        three_LF_MLE_VALUES ["restoreLF"][""];
333        localBranchRes = three_LF_MLES;
334    }
335    writeTheLF (".`thisBranchName`.fit", "three_LF");
336    fprintf (stdout, "\nModel fit:");
337    BIC_scores [thisBranchName] = BIC (localBranchRes[1][0], localBranchRes[1][1]+baseParameters, sample_size);
338    BranchLengthEstimates [thisBranchName] = getBranchLengths ("mixtureTree", 0);
339    fprintf     (csvFilePath, "\nGlobal+",thisBranchName,"\t", localBranchRes[1][0],"\t", localBranchRes[1][1]+baseParameters, "\t", BIC_scores [thisBranchName], "\t", PostOrderAVL2StringDistances(mixTreeAVL,BranchLengthEstimates [thisBranchName]));
340    fprintf (stdout, "\nLocal branch omega distribution: ");
341    reportOmegaDistro(thisBranchName);
342    fprintf (stdout, "Global omega distribution on the rest of the branches: ");
343    reportOmegaDistro ("");
344    stashBranchValues (thisBranchName, "branchValues");
345    pv = 1-CChi2 (2*(localBranchRes[1][0]-res_three_LF_global[1][0]),5);
346    fprintf (stdout, "\nLRT p-value for branch deviation from the global pattern = ", pv, "\n");
347    fprintf     (stdout,"\n", PostOrderAVL2StringDistances(mixTreeAVL,BranchLengthEstimates [thisBranchName]));
348    return 0;
349
350}
351
352//------------------------------------------------------------------------------------------------------------------------
353function writeTheLF (fileNameExt,lfID) {
354    lfOut	= csvFilePath + fileNameExt;//".local.fit";
355    LIKELIHOOD_FUNCTION_OUTPUT = 7;
356    ExecuteCommands ("fprintf (lfOut, CLEAR_FILE, `lfID`);");
357    fprintf (stdout, "[WROTE MODEL FIT TO ", lfOut, "]\n");
358    LIKELIHOOD_FUNCTION_OUTPUT = 2;
359    return 0;
360}
361
362//------------------------------------------------------------------------------------------------------------------------
363function constrainABranch (branch_name, z) {
364    if (z > 0) {
365        ExecuteCommands ("mixtureTree." + branch_name + ".t :< " + z);
366        return 0;
367    } else {
368        if (z < 0) {
369            ExecuteCommands ("mixtureTree." + branch_name + ".t :< 10000");
370            return 0;
371        }
372    }
373    ExecuteCommands ("mixtureTree." + branch_name + ".t = mixtureTreeG." + branch_name + ".t");
374    for (rcc = 1; rcc <= rate_class_count; rcc += 1) {
375        ExecuteCommands ("mixtureTree." + branch_name + ".omega" + rcc + ":= omegaG" + rcc);
376    }
377
378    for (rcc = 1; rcc < rate_class_count; rcc += 1) {
379        ExecuteCommands ("mixtureTree." + branch_name + ".Paux" + rcc + ":= PauxG" + rcc);
380    }
381
382    return 0;
383}
384
385//------------------------------------------------------------------------------------------------------------------------
386function stashBranchValues (branch_name, storage&) {
387
388    locals = {"0" : "t"};
389
390    for (rcc = 1; rcc <= rate_class_count; rcc += 1) {
391        locals + ("omega" + rcc);
392        if (rcc < rate_class_count) {
393            locals + ("Paux" + rcc);
394        }
395    }
396    for (_varID = 0; _varID < Abs (locals); _varID += 1) {
397        storage ["mixtureTree." + branch_name + "." + locals[_varID]] = Eval ("mixtureTree." + branch_name + "." + locals[_varID]);
398    }
399
400    return 0;
401}
402
403//------------------------------------------------------------------------------------------------------------------------
404function unConstrainABranch (branch_name) {
405
406    ExecuteCommands ("mixtureTree." + branch_name + ".t = mixtureTreeG." + branch_name + ".t");
407    for (rcc = 1; rcc <= rate_class_count; rcc += 1) {
408        ExecuteCommands ("mixtureTree." + branch_name + ".omega" + rcc + "= omegaG" + rcc);
409    }
410
411    for (rcc = 1; rcc < rate_class_count; rcc += 1) {
412        ExecuteCommands ("mixtureTree." + branch_name + ".Paux" + rcc + "= PauxG" + rcc);
413    }
414
415    return 0;
416}
417
418//------------------------------------------------------------------------------------------------------------------------
419function BIC (logL, params, sample_size) {
420    BICvalue = -2*logL + 2*params*(sample_size-params-1);
421    fprintf (stdout, " log(L) = ", logL, " with ", params, " parameters, yielding BIC = ", BICvalue, " assuming the sample size of ", sample_size, ".");
422    return BICvalue;
423}
424
425//------------------------------------------------------------------------------------------------------------------------
426function reportOmegaDistro (branchName) {
427    if (Abs(branchName)==0) {
428        for (rcc = 1; rcc <= rate_class_count; rcc += 1) {
429            fprintf (stdout, "\n\t omega = ", Eval("omegaG"+rcc), ", p = ", Eval(global_freqsG[rcc-1]));
430        }
431    } else {
432        for (rcc = 1; rcc < rate_class_count; rcc += 1) {
433            ExecuteCommands ("Paux" + rcc + "=mixtureTree.`branchName`.Paux" +rcc);
434        }
435
436        for (rcc = 1; rcc <= rate_class_count; rcc += 1) {
437            fprintf (stdout, "\n\t omega = ", Eval("mixtureTree.`branchName`.omega"+rcc), ", p = ", Eval(global_freqs[rcc-1]));
438        }
439    }
440    fprintf (stdout, "\n");
441    return 0;
442}
443
444//------------------------------------------------------------------------------------------------------------------------
445
446function getBranchLengthExpression (modelType) {
447    bl_expression = ""; bl_expression * 128;
448
449    if (modelType == 0) {
450        modelPrefix = "MGMatrix";
451        freqMatrix  = global_freqs;
452    } else {
453        modelPrefix = "MGMatrixG";
454        freqMatrix  = global_freqsG;
455    }
456    for (rcc = 1; rcc <= rate_class_count; rcc += 1) {
457        if (rcc > 1) {
458            bl_expression * "+";
459        }
460        ExecuteCommands ("Model _temp = (`modelPrefix`" + rcc + ", codon3x4, 0);");
461        GetString (modelBL, _temp, -1);
462        bl_expression * (freqMatrix[rcc-1]+"*(`modelBL`)");
463    }
464    bl_expression * 0;
465    return bl_expression;
466}
467
468//------------------------------------------------------------------------------------------------------------------------
469
470function getBranchLengths (treeID, modelType) {
471    blByName        = {};
472    bnames_list     = Eval ("BranchName (`treeID`,-1)");
473
474    bl_expression = getBranchLengthExpression(modelType);
475
476    if (modelType == 0) {
477        locals = {"0" : "t"};
478
479        for (rcc = 1; rcc <= rate_class_count; rcc += 1) {
480            locals + ("omega" + rcc);
481            if (rcc < rate_class_count) {
482                locals + ("Paux" + rcc);
483            }
484        }
485
486        for (_bID = 0; _bID < Columns (bnames_list) - 1; _bID += 1) {
487            branch_name = bnames_list[_bID];
488
489             for (_varID = 0; _varID < Abs (locals); _varID += 1) {
490                ExecuteCommands (locals[_varID] + " = Eval (\"`treeID`.`branch_name`." + locals[_varID] + "\")");
491            }
492            blByName [branch_name] = Eval (bl_expression)/3;
493        }
494    } else {
495        for (_bID = 0; _bID < Columns (bnames_list) - 1; _bID += 1) {
496            branch_name = bnames_list[_bID];
497            t = Eval ("`treeID`.`branch_name`.t");
498            blByName [branch_name] = Eval (bl_expression)/3;
499        }
500    }
501
502    return blByName;
503}
504
505
506//------------------------------------------------------------------------------------------------------------------------
507function saveLF (ID)
508{
509	ExecuteCommands ("GetString(_lfInfo,"+ID+",-1)");
510	_stashLF = {};
511	for (_k = 0; _k < Columns (_lfInfo["Global Independent"]); _k+=1)
512	{
513		_stashLF [(_lfInfo["Global Independent"])[_k]] = Eval ((_lfInfo["Global Independent"])[_k]);
514	}
515	for (_k = 0; _k < Columns (_lfInfo["Local Independent"]); _k+=1)
516	{
517		_stashLF [(_lfInfo["Local Independent"])[_k]] = Eval ((_lfInfo["Local Independent"])[_k]);
518	}
519
520	return _stashLF;
521}
522
523//------------------------------------------------------------------------------------------------------------------------
524function restoreLF (key, value)
525{
526	ExecuteCommands (key + " = " + value);
527	return 0;
528}
529
530
531//------------------------------------------------------------------------------------------------------------------------
532function generate_bs_rel_gdd_freqs (numberOfRates, freqs&, model_mixing&, probPrefix, modelPrefix, is_global)
533{
534
535	freqs    	 = {numberOfRates,1};
536	model_mixing	 = ""; model_mixing * 128;
537
538	if (numberOfRates == 1) {
539		freqs[0] = "1";
540	}
541	else {
542	    global_prefix = "";
543
544	    if (is_global)
545	        { global_prefix = "global "; }
546
547        for (mi=1; mi<numberOfRates; mi += 1) {
548            ExecuteCommands (global_prefix+probPrefix+mi+":<1;"+probPrefix+mi+" = 1/" + (numberOfRates-mi+1));
549            ExecuteCommands (global_prefix+probPrefix+mi+":>0;");
550        }
551
552		freqs[0] 	 = ""+probPrefix+"1";
553		for (mi=1; mi<numberOfRates-1; mi+=1) {
554			freqs[mi] = "";
555			for (mi2=1;mi2<=mi;mi2+=1) {
556				freqs[mi] = freqs[mi] + "(1-"+probPrefix+mi2+")*";
557			}
558			freqs[mi] = freqs[mi] + probPrefix+(mi+1);
559		}
560
561		freqs[mi] = "";
562		for (mi2=1;mi2<mi;mi2+=1)
563		{
564			freqs[mi] = freqs[mi] + "(1-"+probPrefix+mi2+")*";
565		}
566		freqs[mi] = freqs[mi] + "(1-"+probPrefix+mi+")";
567	}
568
569	model_mixing * ("Exp(`modelPrefix`1)*"+freqs[0]);
570	for (mi = 1; mi < numberOfRates; mi=mi+1) {
571		model_mixing * ("+Exp(`modelPrefix`"+(mi+1)+")*" + freqs[mi]);
572	}
573	model_mixing * 0;
574	return 0;
575}
576