Skip to content

Commit

Permalink
2.5.60 RC
Browse files Browse the repository at this point in the history
  • Loading branch information
@spond
spond committed Mar 1, 2024
1 parent 82e2945 commit ed17e64
Show file tree
Hide file tree
Showing 13 changed files with 416 additions and 58 deletions.
103 changes: 83 additions & 20 deletions res/TemplateBatchFiles/SelectionAnalyses/BUSTED.bf
View file
Original file line number Diff line number Diff line change
Expand Up @@ -555,8 +555,8 @@ for (busted.partition_index = 0; busted.partition_index < busted.partition_count
busted.global_scaler_list [busted.partition_index] = "busted.bl.scaler_" + busted.partition_index;
parameters.DeclareGlobalWithRanges (busted.global_scaler_list [busted.partition_index], 3, 0, 1000);
busted.initial_ranges [busted.global_scaler_list [busted.partition_index]] = {
terms.lower_bound : 2.5,
terms.upper_bound : 3.5
terms.lower_bound : 2.0,
terms.upper_bound : 4.0
};
}

Expand All @@ -569,7 +569,7 @@ busted.initial_grid = estimators.LHC (busted.initial_ranges,busted.initi
//estimators.CreateInitialGrid (busted.initial_grid, busted.initial_grid.N, busted.initial_grid_presets);

busted.initial_grid = utility.Map (busted.initial_grid, "_v_",
'busted._renormalize (_v_, "busted.distribution", busted.initial.test_mean, busted.error_sink)'
'busted._renormalize_with_weights (_v_, "busted.distribution", busted.initial.test_mean, busted.error_sink)'
);

if (busted.has_background) { //GDD rate category
Expand Down Expand Up @@ -611,24 +611,29 @@ debug.checkpoint = utility.GetEnvVariable ("DEBUG_CHECKPOINT");
if (Type (debug.checkpoint) != "String") {

// constrain nucleotide rate parameters
busted.tmp_fixed = models.FixParameterSetRegExp (terms.nucleotideRatePrefix,busted.test.bsrel_model);
// copy global nucleotide parameter estimates to busted.test.bsrel_model)


busted.tmp_fixed = models.FixParameterSetRegExpFromReference (terms.nucleotideRatePrefix,busted.test.bsrel_model, busted.final_partitioned_mg_results[terms.global]);

busted.grid_search.results = estimators.FitLF (busted.filter_names, busted.trees, busted.model_map, busted.final_partitioned_mg_results, busted.model_object_map,
{
terms.run_options.retain_lf_object: TRUE,
terms.run_options.proportional_branch_length_scaler :
busted.global_scaler_list,

terms.run_options.optimization_settings :
{
"OPTIMIZATION_METHOD" : "nedler-mead",
"MAXIMUM_OPTIMIZATION_ITERATIONS" : 500,
"OPTIMIZATION_PRECISION" : busted.nm.precision
} ,

terms.search_grid : busted.initial_grid,
terms.search_restarts : busted.N.initial_guesses
}
);


busted.grid_search.results = estimators.FitLF (busted.filter_names, busted.trees, busted.model_map, busted.final_partitioned_mg_results, busted.model_object_map, {
"retain-lf-object": TRUE,
terms.run_options.proportional_branch_length_scaler :
busted.global_scaler_list,

terms.run_options.optimization_settings :
{
"OPTIMIZATION_METHOD" : "nedler-mead",
"MAXIMUM_OPTIMIZATION_ITERATIONS" : 500,
"OPTIMIZATION_PRECISION" : busted.nm.precision
} ,

terms.search_grid : busted.initial_grid,
terms.search_restarts : busted.N.initial_guesses
});

parameters.RemoveConstraint (busted.tmp_fixed );

Expand Down Expand Up @@ -1254,6 +1259,64 @@ lfunction busted._renormalize (v, distro, mean, skip_first) {

//------------------------------------------------------------------------------

lfunction busted._renormalize_with_weights (v, distro, mean, skip_first) {

parameters.SetValues (v);
m = parameters.GetStickBreakingDistribution (^distro);
d = Rows (m);

if (skip_first) {
m0 = m[0][0]*m[0][1];
} else {
m0 = 0;
}

over_one = m[d-1][0] * m[d-1][1];

if (over_one >= mean*0.95) {
//console.log ("OVERAGE");
new_weight = mean * Random (0.9, 0.95) / m[d-1][0];
diff = (m[d-1][1] - new_weight)/(d-1);
for (k = (skip_first != 0); k < d-1; k += 1) {
m[k][1] += diff;
}
m[d-1][1] = new_weight;
}

over_one = m[d-1][0] * m[d-1][1];
under_one = (+(m[-1][0] $ m[-1][1])) / (1-m[d-1][1]); // current mean

for (i = (skip_first != 0); i < d-1; i+=1) {
m[i][0] = m[i][0] * mean / under_one;
}

under_one = +(m[-1][0] $ m[-1][1]);

for (i = (skip_first != 0); i < d; i+=1) {
m[i][0] = m[i][0] * mean / under_one;
}

m = m%0;
wts = parameters.SetStickBreakingDistributionWeigths (m[-1][1]);

//console.log (v);

for (i = (skip_first != 0); i < d; i+=1) {
(v[((^distro)["rates"])[i]])[^"terms.fit.MLE"] = m[i][0];
}
for (i = (skip_first != 0); i < d-1; i+=1) {
(v[((^distro)["weights"])[i]])[^"terms.fit.MLE"] = wts[i];
}

//console.log (v);

//assert (0);
return v;

}

//------------------------------------------------------------------------------

lfunction busted.get_multi_hit (model_fit) {
params = selection.io.extract_global_MLE_re (model_fit, '^' + ^'terms.parameters.multiple_hit_rate');
for (k,v; in; selection.io.extract_global_MLE_re (model_fit, '^' + ^'terms.parameters.triple_hit_rate')) {
Expand Down
78 changes: 75 additions & 3 deletions res/TemplateBatchFiles/SelectionAnalyses/RELAX.bf
View file
Original file line number Diff line number Diff line change
Expand Up @@ -519,6 +519,8 @@ if (relax.model_set == "All") { // run all the models
relax.ge_model_map [relax.index] = {"DEFAULT" : "relax.ge"};
}



if (Type (relax.ge_guess) != "Matrix") {


Expand Down Expand Up @@ -677,6 +679,7 @@ if (relax.model_set == "All") { // run all the models
}

} else {
// MINIMAL models branch
relax.general_descriptive.fit = relax.final_partitioned_mg_results;
}

Expand Down Expand Up @@ -833,10 +836,16 @@ if (relax.do_srv) {
if (relax.model_set != "All") {
relax.do_lhc = TRUE;
relax.distribution = models.codon.BS_REL.ExtractMixtureDistribution(relax.model_object_map[relax.reference_model_namespace]);
relax.initial.test_mean = ((selection.io.extract_global_MLE_re (relax.final_partitioned_mg_results, "^" + terms.parameters.omega_ratio + ".+`relax.test_branches_name`.+"))["0"])[terms.fit.MLE];
relax.initial.test_mean = ((selection.io.extract_global_MLE_re (relax.final_partitioned_mg_results, "^" + terms.parameters.omega_ratio + ".+`relax.reference_branches_name`.+"))["0"])[terms.fit.MLE];
relax.init_grid_setup (relax.distribution);
relax.initial_ranges[relax.relaxation_parameter] = {terms.lower_bound : 0.10,
terms.upper_bound : 2.0};


relax.initial_grid = estimators.LHC (relax.initial_ranges,relax.initial_grid.N);
//parameters.SetStickBreakingDistribution (relax.distribution, relax.ge_guess);
relax.initial_grid = utility.Map (relax.initial_grid, "_v_",
'relax._renormalize_with_weights (_v_, "relax.distribution", relax.initial.test_mean)'
);
}

if (relax.has_unclassified) {
Expand Down Expand Up @@ -1001,6 +1010,8 @@ function relax.FitMainTestPair (prompt) {
parameters.DeclareGlobalWithRanges (relax.scaler.id, 1, 0, 1000);
}



relax.general_descriptive.fit = estimators.FitLF (relax.filter_names, relax.trees, relax.model_map,
relax.general_descriptive.fit,
relax.model_object_map,
Expand Down Expand Up @@ -1082,13 +1093,16 @@ function relax.FitMainTestPair (prompt) {


io.ReportProgressMessageMD("RELAX", "alt-2", "* Potential convergence issues due to flat likelihood surfaces; checking to see whether K > 1 or K < 1 is robustly inferred");

if (relax.fitted.K > 1) {
parameters.SetRange (model.generic.GetGlobalParameter (relax.model_object_map ["relax.test"] , terms.relax.k), terms.range01);
} else {
parameters.SetRange (model.generic.GetGlobalParameter (relax.model_object_map ["relax.test"] , terms.relax.k), terms.relax.k_range1);
}

//assert (__SIGTRAP__);

VERBOSITY_LEVEL = 10;

relax.alternative_model.fit.take2 = estimators.FitLF (relax.filter_names, relax.trees, relax.model_map,
relax.alternative_model.fit ,
Expand Down Expand Up @@ -1788,7 +1802,7 @@ function relax.init_grid_setup (omega_distro) {
} else {
relax.initial_ranges [_name_] = {
terms.lower_bound : 1,
terms.upper_bound : 10
terms.upper_bound : 5
};
}
'
Expand Down Expand Up @@ -1846,6 +1860,64 @@ lfunction relax._renormalize (v, distro, mean) {

}

//------------------------------------------------------------------------------

lfunction relax._renormalize_with_weights (v, distro, mean) {

parameters.SetValues (v);
m = parameters.GetStickBreakingDistribution (^distro);
//console.log (v);
//console.log (m);
//console.log (mean);
d = Rows (m);

over_one = m[d-1][0] * m[d-1][1];

if (over_one >= mean*0.95) {
//console.log ("OVERAGE");
new_weight = mean * Random (0.9, 0.95) / m[d-1][0];
diff = (m[d-1][1] - new_weight)/(d-1);
for (k = 0; k < d-1; k += 1) {
m[k][1] += diff;
}
m[d-1][1] = new_weight;
}

over_one = m[d-1][0] * m[d-1][1];
under_one = (+(m[-1][0] $ m[-1][1])) / (1-m[d-1][1]); // current mean

for (i = 0; i < d-1; i+=1) {
m[i][0] = m[i][0] * mean / under_one;
}

under_one = +(m[-1][0] $ m[-1][1]);

for (i = 0; i < d; i+=1) {
m[i][0] = m[i][0] * mean / under_one;
}

m = m%0;
wts = parameters.SetStickBreakingDistributionWeigths (m[-1][1]);



//console.log (v);

for (i = 0; i < d; i+=1) {
(v[((^distro)["rates"])[i]])[^"terms.fit.MLE"] = m[i][0];
}
for (i = 0; i < d-1; i+=1) {
(v[((^distro)["weights"])[i]])[^"terms.fit.MLE"] = wts[i];
}

//console.log (v);

//assert (0);
return v;

}


//------------------------------------------------------------------------------
function relax.optimization_log_file (extension) {
if (relax.OPTIMIZATION_LOGS) {
Expand Down
3 changes: 2 additions & 1 deletion res/TemplateBatchFiles/libv3/all-terms.bf
View file
Original file line number Diff line number Diff line change
Expand Up @@ -396,6 +396,7 @@ namespace terms{
empirical = "empirical";

branch_length_constrain = "branch length constrain";// TODO
branch_length_override = "branch length override";// TODO
get_branch_length = "get-branch-length";
set_branch_length = "set-branch-length";
constrain_branch_length = "constrain-branch-length";
Expand Down Expand Up @@ -509,7 +510,7 @@ namespace terms{
optimization_settings = "optimization-settings";
keep_filters = "use-filters";
optimization_log = "optimization_log";

maintain_branch_lengths_for_grid = "maintain-branch-lengths-for-grid";
}


Expand Down
4 changes: 3 additions & 1 deletion res/TemplateBatchFiles/libv3/models/codon/MG_REV.bf
View file
Original file line number Diff line number Diff line change
Expand Up @@ -244,8 +244,10 @@ function models.codon.MG_REV.set_branch_length(model, value, parameter) {
}
if (utility.Array1D (models.codon.MG_REV.set_branch_length.params)) {
bl_string = Simplify (model[terms.model.branch_length_string],models.codon.MG_REV.set_branch_length.params.subs);

ConvertBranchLength (models.codon.MG_REV.set_branch_length.lp, bl_string , ^(models.codon.MG_REV.set_branch_length.params[0]), 3*value);

ExecuteCommands("FindRoot (models.codon.MG_REV.set_branch_length.lp,(" + bl_string + ")-(" + 3*value + ")," + models.codon.MG_REV.set_branch_length.params[0] + ",0,10000)");
//ExecuteCommands("FindRoot (models.codon.MG_REV.set_branch_length.lp,(" + bl_string + ")-(" + 3*value + ")," + models.codon.MG_REV.set_branch_length.params[0] + ",0,10000)");

Eval (models.codon.MG_REV.set_branch_length.params[0] + "=" + models.codon.MG_REV.set_branch_length.lp);
}
Expand Down
45 changes: 39 additions & 6 deletions res/TemplateBatchFiles/libv3/models/model_functions.bf
View file
Original file line number Diff line number Diff line change
Expand Up @@ -408,7 +408,6 @@ function models.generic.SetBranchLength (model, value, parameter) {
if (Abs((model[terms.parameters])[terms.local]) > 1) {

models.generic.SetBranchLength.bl = Call (model [terms.model.time], model[terms.model.type]);


models.generic.SetBranchLength.bl.p = parameter + "." + models.generic.SetBranchLength.bl;
if (parameters.IsIndependent (models.generic.SetBranchLength.bl.p) == FALSE) {
Expand Down Expand Up @@ -447,33 +446,36 @@ function models.generic.SetBranchLength (model, value, parameter) {
}


if (parameters.IsIndependent (models.generic.SetBranchLength.bl.p)) {
if (parameters.IsIndependent (models.generic.SetBranchLength.bl.p) || model[terms.model.branch_length_override]) {

if (Type (value) == "AssociativeList") {
if (Abs (models.generic.SetBranchLength.expression)) {
ExecuteCommands ("FindRoot (models.generic.SetBranchLength.t,(" + models.generic.SetBranchLength.expression + ")-(" + value[terms.branch_length] + ")," + models.generic.SetBranchLength.bl + ",0,10000)");
ConvertBranchLength (models.generic.SetBranchLength.t, models.generic.SetBranchLength.expression, ^models.generic.SetBranchLength.bl, value[terms.branch_length]);
Eval (models.generic.SetBranchLength.bl.p + ":=(" + value[terms.model.branch_length_scaler] + ")*" + models.generic.SetBranchLength.t);
messages.log ("models.generic.SetBranchLength: " + models.generic.SetBranchLength.bl.p + ":=(" + value[terms.model.branch_length_scaler] + ")*" + models.generic.SetBranchLength.t);

} else {
Eval (models.generic.SetBranchLength.bl.p + ":=(" + value[terms.model.branch_length_scaler] + ")*" + value[terms.branch_length]);
messages.log ("models.generic.SetBranchLength: " + models.generic.SetBranchLength.bl.p + ":=(" + value[terms.model.branch_length_scaler] + ")*" + models.generic.SetBranchLength.t);
//console.log ("models.generic.SetBranchLength: " + models.generic.SetBranchLength.bl.p + ":=(" + value[terms.model.branch_length_scaler] + ")*" + models.generic.SetBranchLength.t);
}

return 1;

} else {

ExecuteCommands ("FindRoot (models.generic.SetBranchLength.t,(" + models.generic.SetBranchLength.expression + ")-(" + value + ")," + models.generic.SetBranchLength.bl + ",0,10000)");
ConvertBranchLength (models.generic.SetBranchLength.t, models.generic.SetBranchLength.expression, ^models.generic.SetBranchLength.bl, value);

//ExecuteCommands ("FindRoot (models.generic.SetBranchLength.t,(" + models.generic.SetBranchLength.expression + ")-(" + value + ")," + models.generic.SetBranchLength.bl + ",0,10000)");
Eval (models.generic.SetBranchLength.bl.p + "=" + models.generic.SetBranchLength.t);
messages.log ("models.generic.SetBranchLength: " + models.generic.SetBranchLength.bl.p + "=" + models.generic.SetBranchLength.t);
}
} else {
messages.log (models.generic.SetBranchLength.bl.p + " was already constrained in models.generic.SetBranchLength");

}
} else {
messages.log ("models.generic.SetBranchLength: missing branch-length-string");
}
}
} else {
messages.log ("models.generic.SetBranchLength: no local model parameters");
}
Expand Down Expand Up @@ -715,3 +717,34 @@ lfunction models.FixParameterSetRegExp (re, model_spec) {
return constraints;
}

/**
* Constrain the set of global model parameters defined by the regexp to own values
* @name parameters.ConstrainParameterSet
* @param {String} re - match this regular expression
* @param {Dict} model_spec - model definition
* @param {Dict} ref - a dictionary of values to use if available
* @returns {Dict} the list of constrained parameters
*/

lfunction models.FixParameterSetRegExpFromReference (re, model_spec, ref) {
constraints = {};
for (tag, id; in; (model_spec[^"terms.parameters"])[^"terms.global"]) {
if (None != regexp.Find (tag, re)) {
if (parameters.IsIndependent (id)) {
mle = ref[tag];
if (Type (mle) == "AssociativeList") {
if (mle / ^"terms.fit.MLE") {
parameters.SetConstraint (id, mle[^"terms.fit.MLE"], "");
constraints + id;
continue;
}
}
parameters.SetConstraint (id, Eval (id), "");
constraints + id;
}
}
}

return constraints;
}

Loading

0 comments on commit ed17e64

Please sign in to comment.