From ed17e64ace4b1efd834bbf71a46fa04068a247f5 Mon Sep 17 00:00:00 2001 From: Sergei Date: Thu, 29 Feb 2024 21:11:50 -0500 Subject: [PATCH] 2.5.60 RC --- .../SelectionAnalyses/BUSTED.bf | 103 ++++++++++++++---- .../SelectionAnalyses/RELAX.bf | 78 ++++++++++++- res/TemplateBatchFiles/libv3/all-terms.bf | 3 +- .../libv3/models/codon/MG_REV.bf | 4 +- .../libv3/models/model_functions.bf | 45 +++++++- .../libv3/models/parameters.bf | 24 ++++ .../libv3/tasks/estimators.bf | 33 ++++-- res/TemplateBatchFiles/libv3/tasks/mpi.bf | 66 +++++++++++ src/core/batchlan.cpp | 36 +++--- src/core/batchlan2.cpp | 6 + src/core/batchlanruntime.cpp | 74 ++++++++++++- src/core/include/batchlan.h | 1 + src/core/include/defines.h | 1 + 13 files changed, 416 insertions(+), 58 deletions(-) diff --git a/res/TemplateBatchFiles/SelectionAnalyses/BUSTED.bf b/res/TemplateBatchFiles/SelectionAnalyses/BUSTED.bf index e711e11c1..c54a0e858 100644 --- a/res/TemplateBatchFiles/SelectionAnalyses/BUSTED.bf +++ b/res/TemplateBatchFiles/SelectionAnalyses/BUSTED.bf @@ -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 }; } @@ -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 @@ -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 ); @@ -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')) { diff --git a/res/TemplateBatchFiles/SelectionAnalyses/RELAX.bf b/res/TemplateBatchFiles/SelectionAnalyses/RELAX.bf index 219128111..4919ed064 100644 --- a/res/TemplateBatchFiles/SelectionAnalyses/RELAX.bf +++ b/res/TemplateBatchFiles/SelectionAnalyses/RELAX.bf @@ -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") { @@ -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; } @@ -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) { @@ -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, @@ -1082,6 +1093,7 @@ 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 { @@ -1089,6 +1101,8 @@ function relax.FitMainTestPair (prompt) { } //assert (__SIGTRAP__); + + VERBOSITY_LEVEL = 10; relax.alternative_model.fit.take2 = estimators.FitLF (relax.filter_names, relax.trees, relax.model_map, relax.alternative_model.fit , @@ -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 }; } ' @@ -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) { diff --git a/res/TemplateBatchFiles/libv3/all-terms.bf b/res/TemplateBatchFiles/libv3/all-terms.bf index 1284b4ba9..5077b4fa9 100644 --- a/res/TemplateBatchFiles/libv3/all-terms.bf +++ b/res/TemplateBatchFiles/libv3/all-terms.bf @@ -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"; @@ -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"; } diff --git a/res/TemplateBatchFiles/libv3/models/codon/MG_REV.bf b/res/TemplateBatchFiles/libv3/models/codon/MG_REV.bf index 1a627e2c1..783bfb5f0 100644 --- a/res/TemplateBatchFiles/libv3/models/codon/MG_REV.bf +++ b/res/TemplateBatchFiles/libv3/models/codon/MG_REV.bf @@ -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); } diff --git a/res/TemplateBatchFiles/libv3/models/model_functions.bf b/res/TemplateBatchFiles/libv3/models/model_functions.bf index 92c89f3bf..c0c9133ca 100644 --- a/res/TemplateBatchFiles/libv3/models/model_functions.bf +++ b/res/TemplateBatchFiles/libv3/models/model_functions.bf @@ -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) { @@ -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"); } @@ -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; +} + diff --git a/res/TemplateBatchFiles/libv3/models/parameters.bf b/res/TemplateBatchFiles/libv3/models/parameters.bf index 16c6e373e..1192dbf1b 100644 --- a/res/TemplateBatchFiles/libv3/models/parameters.bf +++ b/res/TemplateBatchFiles/libv3/models/parameters.bf @@ -727,6 +727,30 @@ lfunction parameters.GetStickBreakingDistributionWeigths (weights) { return w; } +/** + * @name parameters.StStickBreakingDistributionWeigths + * @param {Matrix} 1xN matrix of weights + * @returns {Matrix} 1x(N-1) matrix of weights + */ + +lfunction parameters.SetStickBreakingDistributionWeigths (weights) { + rate_count = utility.Array1D (weights); + + + w = {1, rate_count-1}; + + w[0] = weights[0]; + left_over = (1-w[0]); + + for (i = 1; i < rate_count - 1; i += 1) { + w [i] = weights[i] / left_over; + left_over = left_over * (1-w[i]); + + } + //w[i] = left_over; + return w; +} + /** * @name parameters.helper.stick_breaking diff --git a/res/TemplateBatchFiles/libv3/tasks/estimators.bf b/res/TemplateBatchFiles/libv3/tasks/estimators.bf index f0ff55a11..992ed68a7 100644 --- a/res/TemplateBatchFiles/libv3/tasks/estimators.bf +++ b/res/TemplateBatchFiles/libv3/tasks/estimators.bf @@ -176,7 +176,6 @@ function estimators.SetGlobals2(key2, value) { } - estimators.ApplyExistingEstimates.set_globals[key3] = { terms.id: key3, terms.fit.MLE: value @@ -212,10 +211,11 @@ function estimators.SetGlobals2(key2, value) { * @returns nothing */ function estimators.SetGlobals(key, value) { - /*_did_set = {}; - for (i,v; in; ((value[terms.parameters])[terms.global])) { - _did_set [v] = 0; - }*/ + //_did_set = {}; + //console.log (((value[terms.parameters])[terms.global])); + //for (i,v; in; ((value[terms.parameters])[terms.global])) { + // _did_set [v] = 0; + //} ((value[terms.parameters])[terms.global])["estimators.SetGlobals2"][""]; @@ -496,9 +496,12 @@ function estimators.ApplyExistingEstimatesToTree (_tree_name, model_descriptions estimators.ApplyExistingEstimatesToTree.constraint_count = 0; - ExecuteCommands("GetInformation (estimators.ApplyExistingEstimatesToTree.map, `_tree_name`);"); estimators.ApplyExistingEstimatesToTree.branch_names = Rows(estimators.ApplyExistingEstimatesToTree.map); + + for (i,v; in; model_descriptions) { + parameters.SetCategoryVariables (v); + } for (estimators.ApplyExistingEstimatesToTree.b = 0; estimators.ApplyExistingEstimatesToTree.b < Abs(estimators.ApplyExistingEstimatesToTree.map); estimators.ApplyExistingEstimatesToTree.b += 1) { _branch_name = estimators.ApplyExistingEstimatesToTree.branch_names[estimators.ApplyExistingEstimatesToTree.b]; @@ -565,6 +568,8 @@ function estimators.ApplyExistingEstimates(likelihood_function_id, model_descrip estimators.ApplyExistingEstimates.df_correction = 0; // copy global variables first + + estimators.ApplyExistingEstimates.results[terms.global] = {}; model_descriptions["estimators.SetCategory"][""]; // the above line traverses all model descriptions and sets @@ -591,6 +596,8 @@ function estimators.ApplyExistingEstimates(likelihood_function_id, model_descrip _application_type = None; + //console.log (branch_length_conditions); + if (Type (branch_length_conditions) == "AssociativeList") { if (Abs(branch_length_conditions) > estimators.ApplyExistingEstimates.i) { _application_type = branch_length_conditions[estimators.ApplyExistingEstimates.i]; @@ -762,19 +769,24 @@ lfunction estimators.FitLF(data_filter, tree, model_map, initial_values, model_o if (Type(initial_values) == "AssociativeList") { utility.ToggleEnvVariable("USE_LAST_RESULTS", 1); - //console.log (initial_values); df = estimators.ApplyExistingEstimates("`&likelihoodFunction`", model_objects, initial_values, run_options[utility.getGlobalValue("terms.run_options.proportional_branch_length_scaler")]); - } - + } + can_do_restarts = null; if (utility.Has (run_options, utility.getGlobalValue("terms.search_grid"),"AssociativeList")) { - grid_results = mpi.ComputeOnGrid (&likelihoodFunction, run_options [utility.getGlobalValue("terms.search_grid")], "mpi.ComputeOnGrid.SimpleEvaluator", "mpi.ComputeOnGrid.ResultHandler"); + grid_results = mpi.ComputeOnGridSetValues (&likelihoodFunction, run_options [utility.getGlobalValue("terms.search_grid")], { + 0: model_objects, + 1: initial_values, + 2: run_options[utility.getGlobalValue("terms.run_options.proportional_branch_length_scaler")] + }, "mpi.ComputeOnGrid.SimpleEvaluatorWithValues", "mpi.ComputeOnGrid.ResultHandler"); + if (utility.Has (run_options, utility.getGlobalValue("terms.search_restarts"),"Number")) { restarts = run_options[utility.getGlobalValue("terms.search_restarts")]; if (restarts > 1) { grid_results = utility.DictToSortedArray (grid_results); + //console.log (grid_results); can_do_restarts = {}; for (i = 1; i <= restarts; i += 1) { can_do_restarts + (run_options [utility.getGlobalValue("terms.search_grid")])[grid_results[Rows(grid_results)-i][1]]; @@ -787,6 +799,7 @@ lfunction estimators.FitLF(data_filter, tree, model_map, initial_values, model_o } //console.log (best_value); //console.log ((run_options [utility.getGlobalValue("terms.search_grid")])[best_value["key"]]); + //console.log (can_do_restarts); //assert (0); } diff --git a/res/TemplateBatchFiles/libv3/tasks/mpi.bf b/res/TemplateBatchFiles/libv3/tasks/mpi.bf index c1f5eb384..a94a61d01 100644 --- a/res/TemplateBatchFiles/libv3/tasks/mpi.bf +++ b/res/TemplateBatchFiles/libv3/tasks/mpi.bf @@ -298,6 +298,71 @@ namespace mpi { } } + //------------------------------------------------------------------------------ + + lfunction ComputeOnGridSetValues (lf_id, grid, values, handler, callback) { + + jobs = mpi.PartitionIntoBlocks(grid); + + scores = {}; + + + queue = mpi.CreateQueue ({^"terms.mpi.LikelihoodFunctions": {{lf_id}}, + ^"terms.mpi.Headers" : utility.GetListOfLoadedModules ("libv3/")}); + + io.ReportProgressBar("", "Computing LF on a grid"); + for (i = 1; i < Abs (jobs); i += 1) { + //io.ReportProgressBar("", "Computing LF on a grid " + i + "/" + Abs (jobs)); + mpi.QueueJob (queue, handler, {"0" : lf_id, + "1" : jobs [i], + "2" : values, + "3" : &scores}, callback); + } + + + Call (callback, -1, Call (handler, lf_id, jobs[0], values, &scores), {"0" : lf_id, "1" : jobs [0], "3" : values, "2" : &scores}); + mpi.QueueComplete (queue); + + io.ClearProgressBar(); + + return scores; + } + + + //------------------------------------------------------------------------------ + + lfunction ComputeOnGrid.SimpleEvaluatorWithValues (lf_id, tasks, values, scores) { + + + //#profile START; + + LFCompute (^lf_id, LF_START_COMPUTE); + + results = {}; + task_ids = utility.Keys (tasks); + task_count = Abs (tasks); + for (i, v; in; values[0]) { + v[^"terms.model.branch_length_override"] = TRUE; + } + for (i = 0; i < task_count; i+=1) { + + parameters.SetValues (tasks[task_ids[i]]); + estimators.ApplyExistingEstimates (lf_id, values[0], values[1], values[2]); + //fprintf (stdout, ^lf_id); + LFCompute (^lf_id, ll); + results [task_ids[i]] = ll; + io.ReportProgressBar("", "Grid result " + i + " = " + ll + " (best = " + Max (results, 0)[^"terms.data.value"] + ")"); + + } + LFCompute (^lf_id, LF_DONE_COMPUTE); + for (i, v; in; values[0]) { + v - ^"terms.model.branch_length_override"; + } + + //#profile _hyphy_profile_dump; + //utility.FinishAndPrintProfile (_hyphy_profile_dump); + return results; + } //------------------------------------------------------------------------------ @@ -338,6 +403,7 @@ namespace mpi { //------------------------------------------------------------------------------ lfunction ComputeOnGrid.SimpleEvaluator (lf_id, tasks, scores) { + LFCompute (^lf_id, LF_START_COMPUTE); results = {}; diff --git a/src/core/batchlan.cpp b/src/core/batchlan.cpp index 2a114fb91..f34141369 100644 --- a/src/core/batchlan.cpp +++ b/src/core/batchlan.cpp @@ -2001,6 +2001,7 @@ bool _ExecutionList::BuildList (_StringBuffer& s, _SimpleList* bc, bool case HY_HBL_COMMAND_CONSTRUCT_CATEGORY_MATRIX: case HY_HBL_COMMAND_KEYWORD_ARGUMENT: case HY_HBL_COMMAND_DO_SQL: + case HY_HBL_COMMAND_CONVERT_BRANCH_LENGTH: { _ElementaryCommand::ExtractValidateAddHBLCommand (currentLine, prefixTreeCode, pieces, commandExtraInfo, *this); break; @@ -2168,26 +2169,25 @@ _ElementaryCommand::_ElementaryCommand (_String& s) { //____________________________________________________________________________________ _ElementaryCommand::~_ElementaryCommand (void) { - if (CanFreeMe()) { - if (code==4) { // IF - if (simpleParameters.lLength>2) { - _Formula* f = (_Formula*)simpleParameters(2); + + auto clear_formulas = [this] (long start_at) -> void { + for (long i = start_at; i 2) { delete (AVLListXLIterator*)simpleParameters.get(2); @@ -2341,6 +2341,7 @@ BaseRef _ElementaryCommand::toStr (unsigned long) { case HY_HBL_COMMAND_SELECT_TEMPLATE_MODEL: case HY_HBL_COMMAND_KEYWORD_ARGUMENT: case HY_HBL_COMMAND_GET_INFORMATION: + case HY_HBL_COMMAND_CONVERT_BRANCH_LENGTH: case HY_HBL_COMMAND_SIMULATE_DATA_SET: { (*string_form) << procedure (code); } @@ -3713,6 +3714,9 @@ bool _ElementaryCommand::Execute (_ExecutionList& chain) { return HandleFindRootOrIntegrate(chain, code == HY_HBL_COMMAND_INTEGRATE); break; + case HY_HBL_COMMAND_CONVERT_BRANCH_LENGTH: + return HandleConvertBranchLength(chain); + break; case HY_HBL_COMMAND_MPI_SEND: return HandleMPISend (chain); diff --git a/src/core/batchlan2.cpp b/src/core/batchlan2.cpp index 321d6d027..b76311970 100644 --- a/src/core/batchlan2.cpp +++ b/src/core/batchlan2.cpp @@ -175,6 +175,8 @@ void _HBL_Init_Const_Arrays (void) _HY_ValidHBLExpressions.Insert ("BGM ", HY_HBL_COMMAND_BGM); _HY_ValidHBLExpressions.Insert ("SimulateDataSet", HY_HBL_COMMAND_SIMULATE_DATA_SET); _HY_ValidHBLExpressions.Insert ("KeywordArgument", HY_HBL_COMMAND_KEYWORD_ARGUMENT); + _HY_ValidHBLExpressions.Insert ("ConvertBranchLength(", HY_HBL_COMMAND_CONVERT_BRANCH_LENGTH); + /* const long cut, const long conditions, const char sep, const bool doTrim, const bool isAssignment, const bool needsVerb, length options */ @@ -462,6 +464,10 @@ const long cut, const long conditions, const char sep, const bool doTrim, const // matrix global arrays + _HY_HBLCommandHelper.Insert ((BaseRef)HY_HBL_COMMAND_CONVERT_BRANCH_LENGTH, + (long)_hyInitCommandExtras (_HY_ValidHBLExpressions.Insert ("ConvertBranchLength(", HY_HBL_COMMAND_CONVERT_BRANCH_LENGTH,false), + 4, + "ConvertBranchLength(, , , )",',')); _HY_MatrixRandomValidPDFs.Insert ("Dirichlet", _HY_MATRIX_RANDOM_DIRICHLET, "Gaussian", _HY_MATRIX_RANDOM_GAUSSIAN, diff --git a/src/core/batchlanruntime.cpp b/src/core/batchlanruntime.cpp index 92e29be79..50e5591ea 100644 --- a/src/core/batchlanruntime.cpp +++ b/src/core/batchlanruntime.cpp @@ -314,6 +314,79 @@ bool _ElementaryCommand::HandleDifferentiate(_ExecutionList& current_progra //____________________________________________________________________________________ +bool _ElementaryCommand::HandleConvertBranchLength(_ExecutionList& current_program){ + + _Variable * receptacle = nil; + current_program.advance (); + + auto cleanup_cache = [this] (long start_at) -> void { + for (long k = simpleParameters.countitems() - 1; k >= start_at; k --) { + _Formula* f = (_Formula*)simpleParameters.get (k); + if (f) { + delete f; + } + simpleParameters.Delete(k); + } + }; + + try { + receptacle = _ValidateStorageVariable (current_program); + _List dynamic_manager; + _FString* expression = (_FString*)_ProcessAnArgumentByType(*GetIthParameter(1), STRING, current_program, &dynamic_manager); + _Formula * lhs_expression = nil, + * lhs_derivative = nil; + + if (parameters.countitems () > 4) { + _String * cached_fla = GetIthParameter(4); + if (!cached_fla->Equal(expression->get_str()) || simpleParameters.empty()) { + cleanup_cache (0); + parameters.Delete (4); + } else { + lhs_expression = (_Formula*)simpleParameters.get (0); + if (simpleParameters.countitems() > 1) { + lhs_derivative = (_Formula*)simpleParameters.get (1); + } + } + } + + if (parameters.countitems () <= 4) { + parameters < expression->get_str_ref(); + } + + if (!lhs_expression) { + lhs_expression = new _Formula; + _CheckExpressionForCorrectness (*lhs_expression, expression->get_str(), current_program); + simpleParameters << (long)lhs_expression; + } + + + _Variable * target_variable = _CheckForExistingVariableByType (*GetIthParameter(2),current_program,NUMBER); + + if (!lhs_expression->DependsOnVariable(target_variable->get_index())) { + throw (expression->get_str().Enquote() & " does not depend on the variable " & target_variable->GetName()->Enquote()); + } + + hyFloat target_value = _ProcessNumericArgumentWithExceptions(*GetIthParameter(3),current_program.nameSpacePrefix); + + if (!lhs_derivative && simpleParameters.countitems() < 2) { + lhs_derivative = lhs_expression->Differentiate (*target_variable->GetName(),false); + simpleParameters << (long)lhs_derivative; + } + + if (lhs_derivative) { + receptacle->SetValue (new _Constant (lhs_expression->Newton (*lhs_derivative,target_variable, target_value, 0.0, 10000.)),false,true, NULL); + } else { + receptacle->SetValue (new _Constant (lhs_expression->Brent (target_variable, 0.0, 10000.,1.e-7, nil, target_value)), false,true, NULL); + } + + } catch (const _String& error) { + return _DefaultExceptionHandler (receptacle, error, current_program); + } + return true; +} + +//____________________________________________________________________________________ + bool _ElementaryCommand::HandleFindRootOrIntegrate (_ExecutionList& currentProgram, bool do_integrate){ _Variable * receptacle = nil; @@ -331,7 +404,6 @@ bool _ElementaryCommand::HandleFindRootOrIntegrate (_ExecutionList& current if (!parsed_expression.DependsOnVariable(target_variable->get_index()) && !do_integrate) { throw (expression.Enquote() & " does not depend on the variable " & target_variable->GetName()->Enquote()); } - hyFloat lb = _ProcessNumericArgumentWithExceptions (*GetIthParameter(3),currentProgram.nameSpacePrefix), ub = _ProcessNumericArgumentWithExceptions (*GetIthParameter(4),currentProgram.nameSpacePrefix); diff --git a/src/core/include/batchlan.h b/src/core/include/batchlan.h index 4ed0b20e3..e14d456f0 100644 --- a/src/core/include/batchlan.h +++ b/src/core/include/batchlan.h @@ -282,6 +282,7 @@ class _ElementaryCommand: public _String // string contains the literal for th bool HandleExport (_ExecutionList&); bool HandleDifferentiate (_ExecutionList&); bool HandleFindRootOrIntegrate (_ExecutionList&, bool do_integrate = false); + bool HandleConvertBranchLength (_ExecutionList&); bool HandleMPISend (_ExecutionList&); bool HandleMPIReceive (_ExecutionList&); bool HandleExecuteCommandsCases (_ExecutionList&, bool do_load_from_file = false, bool do_load_library = false); diff --git a/src/core/include/defines.h b/src/core/include/defines.h index 767f81024..18fbe3225 100644 --- a/src/core/include/defines.h +++ b/src/core/include/defines.h @@ -319,6 +319,7 @@ SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. #define HY_HBL_COMMAND_KEYWORD_ARGUMENT 567L #define HY_HBL_COMMAND_INIT_ITERATOR 568L #define HY_HBL_COMMAND_ADVANCE_ITERATOR 569L +#define HY_HBL_COMMAND_CONVERT_BRANCH_LENGTH 570L //! HyPhy standard directory locations