diff --git a/CMakeLists.txt b/CMakeLists.txt index 3812aad9a..3ecbba0a1 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -15,6 +15,7 @@ option(NOAVX OFF) option(NOSSE4 OFF) option(NONEON OFF) option(NOZLIB OFF) +option(NOBLAS OFF) if(CMAKE_SYSTEM_NAME STREQUAL "Emscripten") set(NOAVX ON) @@ -427,7 +428,7 @@ include_directories( contrib/SQLite-3.8.2 # SQLite ) -if (${BLAS_FOUND}) +if (${BLAS_FOUND} AND NOT NOBLAS) add_definitions (-D_SLKP_USE_APPLE_BLAS) set(DEFAULT_LIBRARIES "${DEFAULT_LIBRARIES}") endif() @@ -451,7 +452,7 @@ else(${OPENMP_FOUND}) target_link_libraries(hyphy PRIVATE ${DEFAULT_LIBRARIES}) endif(${OPENMP_FOUND}) -if(${BLAS_FOUND}) +if(${BLAS_FOUND} AND NOT NOBLAS) target_link_libraries(hyphy PRIVATE ${BLAS_LIBRARIES}) endif () @@ -490,7 +491,7 @@ if(${MPI_FOUND}) target_link_libraries(HYPHYMPI ${DEFAULT_LIBRARIES} ${MPI_LIBRARIES}) endif(${OPENMP_FOUND}) - if(${BLAS_FOUND}) + if(${BLAS_FOUND} AND NOT NOBLAS) target_link_libraries(HYPHYMPI ${BLAS_LIBRARIES}) endif () @@ -598,6 +599,9 @@ target_link_libraries(HYPHYDEBUG ${DEFAULT_LIBRARIES}) add_custom_target(DEBUG DEPENDS HYPHYDEBUG) +if(${BLAS_FOUND}) + target_link_libraries(HYPHYDEBUG ${BLAS_LIBRARIES}) +endif () set_target_properties( HYPHYDEBUG diff --git a/res/TemplateBatchFiles/MSS-joint-fitter.bf b/res/TemplateBatchFiles/MSS-joint-fitter.bf index 2c627a644..596dc7837 100644 --- a/res/TemplateBatchFiles/MSS-joint-fitter.bf +++ b/res/TemplateBatchFiles/MSS-joint-fitter.bf @@ -53,6 +53,13 @@ io.ReportProgressMessageMD("mss", "data" , "* Loaded a list with **" + mss_selec KeywordArgument ("code", "Which genetic code should be used", "Universal"); mss.genetic_code = alignments.LoadGeneticCode (None); +KeywordArgument ("omega", "How should alignment-level omega be treated?", "Fix"); + +mss.omega_option = io.SelectAnOption ({ + {"Fix", "[Default] Fix omega estimates at those obtained with the standard MG94xREV model"} + {"Fit", "Fit omega (per alignment) together with the MSS model"} + }, "How should alignment-level omega be treated?"); + mss.file_records = {}; mss.file_info = {}; @@ -160,10 +167,18 @@ mss_selector.header = { {"TreeLength"}, {"Log(L)"} }; + +KeywordArgument ("model", "Substitution model to use","SynREV"); +mss.model_option = io.SelectAnOption ({ + {"SynREV", "[Default] SynREV model (one rate per aa)"} + {"SynREVCodon", "SynREVCodon (one rate per codon pair)"} + }, "Which model?"); + + ExecuteCommands ( "mss.codon_classes = model.codon.MSS.prompt_and_define (terms.global, mss.genetic_code[terms.code])", - {"--mss-type" : "SynREV"} + {"--mss-type" : mss.model_option} ); io.ReportProgressMessageMD("mss", "fit0" , "Individual file statistics and simple model fits\n"); @@ -219,6 +234,7 @@ io.ReportProgressMessageMD("mss", "fit0done", "**LogL = " + mss.baselineLL + "** mss.json ["baseline"] = {terms.json.log_likelihood : mss.baselineLL, terms.json.AICc : mss.baseline_AIC}; function mss.MSS_generator (type) { + mss.codon_classes [utility.getGlobalValue("terms.model.MSS.normalize")] = TRUE; model = Call ("models.codon.MSS.ModelDescription",type, mss.genetic_code[terms.code], mss.codon_classes); return model; } @@ -234,10 +250,13 @@ mss.lf_id = "MSS.COMPOSITE.LF"; mss.model_map_overall = {}; +//#profile START; +mss.constraints = {}; for (mss.counter, mss_selector.path; in; mss_selector.file_list) { mss.prefix = mss.file_prefix [mss_selector.path]; + console.log (">" + mss.counter + " / " + mss_selector.path); mss.model_name = "`mss.prefix`.model"; mss.tree_name = "`mss.prefix`.tree"; @@ -271,14 +290,17 @@ for (mss.counter, mss_selector.path; in; mss_selector.file_list) { (initial_values [terms.branch_length])[0], terms.model.branch_length_constrain, TRUE); - - models.FixParameterSetRegExp (terms.parameters.omega_ratio, ^(mss.model_name)); + + if (mss.omega_option == "Fix") { + models.FixParameterSetRegExp (terms.parameters.omega_ratio, ^(mss.model_name)); + mss.constraint_count += 1; + } models.FixParameterSetRegExp (terms.nucleotideRatePrefix, ^(mss.model_name)); if (mss.counter == 0) { mss.reference_set = ^(mss.model_name); - mss.mss_rate_list = model.GetParameters_RegExp( ^(mss.model_name),"^" + terms.parameters.synonymous_rate + ""); + mss.mss_rate_list = model.GetParameters_RegExp( ^(mss.model_name),"^" + terms.MeanScaler("")); mss.model_dimension = utility.Array1D (mss.mss_rate_list); mss.scaler_prefix = 'mss.scaler_parameter_'; mss.scaler_mapping = {}; @@ -290,15 +312,18 @@ for (mss.counter, mss_selector.path; in; mss_selector.file_list) { for (mss.i2 = 0; mss.i2 < mss.model_dimension; mss.i2 += 1) { parameters.DeclareGlobal (mss.scaler_prefix + mss.i2, None); } - mss.json ["mapping"] = mss.scaler_index; } else { //utility.getGlobalValue ("terms.parameters.synonymous_rate"); - models.BindGlobalParameters ({"0" : mss.reference_set , "1" : ^(mss.model_name)}, terms.parameters.synonymous_rate + terms.model.MSS.syn_rate_between); - models.BindGlobalParameters ({"0" : mss.reference_set , "1" : ^(mss.model_name)}, terms.parameters.synonymous_rate + terms.model.MSS.syn_rate_within); - } + mss.constraints * models.BindGlobalParametersDeferred ({"0" : mss.reference_set , "1" : ^(mss.model_name)}, "^" + terms.MeanScaler("")); + } } +parameters.BatchApplyConstraints (mss.constraints); + +//#profile _hyphy_profile_dump; +//utility.FinishAndPrintProfile (_hyphy_profile_dump); +utility.SetEnvVariable ("AUTO_PARALLELIZE_OPTIMIZE", 1); utility.ExecuteInGlobalNamespace ("LikelihoodFunction `mss.lf_id` = (`&mss.lf_components`)"); @@ -317,6 +342,17 @@ io.ReportProgressMessageMD("mss", "fitAlldone", "**LogL = " + res[1][0] + "**" mss.json["joint-model"] = estimators.ExtractMLEsOptions (mss.lf_id, mss.model_map_overall, {terms.globals_only : TRUE}); //estimators.RemoveBranchLengthConstraints (mss.json["joint-model"]); +function pfilter (key, value) { + if (key[0] != "[" || (key $ "ratio")[0] >= 0) { + mss.filtered [key] = value; + } +} + +mss.filtered = {}; +((mss.json["joint-model"])[terms.global])["pfilter"][""]; +(mss.json ["joint-model"])[terms.global] = mss.filtered; + + io.SpoolJSON (mss.json, mss.json[terms.data.file]); KeywordArgument ("save-fit", "Write the resulting model fit file to this (large!) file", "/dev/null"); diff --git a/res/TemplateBatchFiles/SelectionAnalyses/BUSTED.bf b/res/TemplateBatchFiles/SelectionAnalyses/BUSTED.bf index c54a0e858..e1f513cde 100644 --- a/res/TemplateBatchFiles/SelectionAnalyses/BUSTED.bf +++ b/res/TemplateBatchFiles/SelectionAnalyses/BUSTED.bf @@ -48,7 +48,6 @@ Version 4.5 adds an 'error absorption' component [experimental] }; -console.log (busted.analysis_description); io.DisplayAnalysisBanner (busted.analysis_description); busted.FG = "Test"; @@ -536,8 +535,7 @@ if (busted.do_srv) { //PARAMETER_GROUPING + busted.srv_distribution["weights"]; PARAMETER_GROUPING + utility.Concat (busted.srv_distribution["rates"],busted.srv_distribution["weights"]); - busted.init_grid_setup (busted.srv_distribution, FALSE); - + busted.init_grid_setup (busted.srv_distribution, FALSE); } if (busted.do_srv_hmm) { @@ -553,33 +551,29 @@ busted.global_scaler_list = {}; for (busted.partition_index = 0; busted.partition_index < busted.partition_count; busted.partition_index += 1) { 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.0, - terms.upper_bound : 4.0 - }; + parameters.DeclareGlobalWithRanges (busted.global_scaler_list [busted.partition_index], 1, 0, 1000); + // busted.initial_ranges [busted.global_scaler_list [busted.partition_index]] = { + // terms.lower_bound : 1.0, + // terms.upper_bound : 3.0 + // }; } - - busted.initial.test_mean = ((selection.io.extract_global_MLE_re (busted.final_partitioned_mg_results, "^" + terms.parameters.omega_ratio + ".+test.+"))["0"])[terms.fit.MLE]; busted.initial_grid = estimators.LHC (busted.initial_ranges,busted.initial_grid.N); -//estimators.CreateInitialGrid (busted.initial_grid, busted.initial_grid.N, busted.initial_grid_presets); - busted.initial_grid = utility.Map (busted.initial_grid, "_v_", 'busted._renormalize_with_weights (_v_, "busted.distribution", busted.initial.test_mean, busted.error_sink)' ); -if (busted.has_background) { //GDD rate category +if (busted.has_background) { //has background busted.initial.background_mean = ((selection.io.extract_global_MLE_re (busted.final_partitioned_mg_results, "^" + terms.parameters.omega_ratio + ".+background.+"))["0"])[terms.fit.MLE]; + busted.initial_grid = utility.Map (busted.initial_grid, "_v_", - 'busted._renormalize (_v_, "busted.background_distribution", busted.initial.background_mean, busted.error_sink)' + 'busted._renormalize_with_weights (_v_, "busted.background_distribution", busted.initial.background_mean, busted.error_sink)' ); } - busted.model_map = {}; for (busted.partition_index = 0; busted.partition_index < busted.partition_count; busted.partition_index += 1) { @@ -604,7 +598,7 @@ io.ReportProgressMessageMD ("BUSTED", "main", "Performing the full (dN/dS > 1 al */ -busted.nm.precision = -0.00025*busted.final_partitioned_mg_results[terms.fit.log_likelihood]; +busted.nm.precision = Max (-0.00001*busted.final_partitioned_mg_results[terms.fit.log_likelihood],0.5); debug.checkpoint = utility.GetEnvVariable ("DEBUG_CHECKPOINT"); @@ -615,7 +609,8 @@ if (Type (debug.checkpoint) != "String") { 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, @@ -634,8 +629,12 @@ if (Type (debug.checkpoint) != "String") { } ); + + parameters.RemoveConstraint (busted.tmp_fixed ); + //console.log (busted.tmp_fixed); + PARAMETER_GROUPING + Columns (busted.tmp_fixed); //PRODUCE_OPTIMIZATION_LOG = 1; @@ -643,7 +642,7 @@ if (Type (debug.checkpoint) != "String") { "retain-lf-object": TRUE, terms.run_options.optimization_settings : { - "OPTIMIZATION_METHOD" : "coordinate-wise", + "OPTIMIZATION_METHOD" : "hybrid", //"OPTIMIZATION_PRECISION" : 1. } @@ -1260,46 +1259,68 @@ 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); - + + mean = Max (mean, 1e-3); + 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); + diff = (m[d-1][1] - new_weight)/(d-1-(skip_first != 0)); 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 + under_one = (+(m[-1][0] $ m[-1][1]) - m0) / (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]); + + if (skip_first) { + m_rest = m [{{1,0}}][{{d-1,1}}]; + under_one = +(m_rest[-1][0] $ m_rest[-1][1]); + } else { + 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]); + + + if (skip_first) { + m_rest = m [{{1,0}}][{{d-1,1}}]; + m_rest = m_rest % 0; + for (i = 1; i < d; i+=1) { + m[i][0] = m_rest[i-1][0]; + m[i][1] = m_rest[i-1][1]; + } + } else { + m = m%0; + } - //console.log (v); + + wts = parameters.SetStickBreakingDistributionWeigths (m[-1][1]); + for (i = (skip_first != 0); i < d; i+=1) { (v[((^distro)["rates"])[i]])[^"terms.fit.MLE"] = m[i][0]; @@ -1307,10 +1328,6 @@ lfunction busted._renormalize_with_weights (v, distro, mean, skip_first) { 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; } @@ -1383,7 +1400,7 @@ function busted.init_grid_setup (omega_distro, error_sink) { busted.initial_grid [_name_] = {{100,500,1000,5000}}; busted.initial_ranges [_name_] = { terms.lower_bound : 100, - terms.upper_bound : 10000 + terms.upper_bound : 1000 }; } else { busted.initial_ranges [_name_] = { @@ -1416,12 +1433,12 @@ function busted.init_grid_setup (omega_distro, error_sink) { if (error_sink && _index_ == 0) { busted.initial_grid [_name_] = { { - 0, 0.001, 0.005, 0.1 + 0, 0.001, 0.0025, 0.025 } }; busted.initial_ranges [_name_] = { terms.lower_bound : 0, - terms.upper_bound : 0.01, + terms.upper_bound : 0.005, }; } else { busted.initial_ranges [_name_] = { diff --git a/res/TemplateBatchFiles/SelectionAnalyses/FEL.bf b/res/TemplateBatchFiles/SelectionAnalyses/FEL.bf index e3ed69629..d82f5ab18 100644 --- a/res/TemplateBatchFiles/SelectionAnalyses/FEL.bf +++ b/res/TemplateBatchFiles/SelectionAnalyses/FEL.bf @@ -270,7 +270,7 @@ if (^"fel.ci") { } } else { - if (^"fel.resample" > 0 ) { + if (^"fel.resample" == 0 ) { if (fel.multi_hit == "Double") { fel.table_headers = { {"alpha", "Synonymous substitution rate at a site"} @@ -980,7 +980,7 @@ lfunction fel.store_results (node, result, arguments) { if (^"fel.ci") { - if (^"fel.resample") { + if (^"fel.resample" > 0) { result_row = { 0 : 0, // alpha 1 : 0 , // beta 2: 0, // alpha==beta @@ -1014,6 +1014,7 @@ lfunction fel.store_results (node, result, arguments) { 3:0, // LRT 4:1, // p-value, 5:0, // total branch length of tested branches + 6:0 // asymptotic p-value } ; } else { result_row = { 0:0, // alpha @@ -1040,15 +1041,17 @@ lfunction fel.store_results (node, result, arguments) { has_lrt = result / ^"terms.simulated"; + + if (^"fel.resample") { - result_row [utility.Array1D (result_row) - 1] = lrt [utility.getGlobalValue("terms.p_value")]; + result_row [N_col - 1] = lrt [utility.getGlobalValue("terms.p_value")]; } - + if (has_lrt) { pv = +((result[^"terms.simulated"])["_MATRIX_ELEMENT_VALUE_>=" + lrt [utility.getGlobalValue("terms.LRT")]]); lrt [utility.getGlobalValue("terms.p_value")] = (pv+1)/(1+^"fel.resample"); } - + result_row [0] = estimators.GetGlobalMLE (result[utility.getGlobalValue("terms.alternative")], ^"fel.site_alpha"); result_row [1] = estimators.GetGlobalMLE (result[utility.getGlobalValue("terms.alternative")], ^"fel.site_beta"); result_row [2] = estimators.GetGlobalMLE (result[utility.getGlobalValue("terms.Null")], ^"fel.site_beta"); diff --git a/res/TemplateBatchFiles/SelectionAnalyses/RELAX.bf b/res/TemplateBatchFiles/SelectionAnalyses/RELAX.bf index 10015b390..166eee840 100644 --- a/res/TemplateBatchFiles/SelectionAnalyses/RELAX.bf +++ b/res/TemplateBatchFiles/SelectionAnalyses/RELAX.bf @@ -836,19 +836,15 @@ 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.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); - relax.initial_grid = utility.Map (relax.initial_grid, "_v_", - 'relax._renormalize_with_weights (_v_, "relax.distribution", relax.initial.test_mean)' - ); } if (relax.has_unclassified) { + relax.unclassified.bsrel_model = model.generic.DefineMixtureModel(relax.model_generator, "relax.unclassified", { "0": parameters.Quote(terms.global), @@ -866,19 +862,39 @@ if (relax.has_unclassified) { relax.distribution_uc = models.codon.BS_REL.ExtractMixtureDistribution(relax.unclassified.bsrel_model); relax.init_grid_setup (relax.distribution_uc); } + parameters.SetRange (model.generic.GetGlobalParameter (relax.unclassified.bsrel_model , terms.AddCategory (terms.parameters.omega_ratio,relax.rate_classes)), terms.range_gte1); relax.model_object_map ["relax.unclassified"] = relax.unclassified.bsrel_model; - - for (relax.index, relax.junk ; in; relax.filter_names) { + for (relax.index, relax.junk ; in; relax.filter_names) { (relax.model_map[relax.index]) ["relax.unclassified"] = utility.Filter (relax.selected_branches[relax.index], '_value_', '_value_ == relax.unclassified_branches_name'); } models.BindGlobalParameters ({"0" : relax.model_object_map[relax.reference_model_namespace], "1" : relax.unclassified.bsrel_model}, terms.nucleotideRate("[ACGT]","[ACGT]")); } +if (relax.do_lhc) { + relax.initial_grid = estimators.LHC (relax.initial_ranges,relax.initial_grid.N); + 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.initial_grid = utility.Map (relax.initial_grid, "_v_", + 'relax._renormalize_with_weights (_v_, "relax.distribution", relax.initial.test_mean)' + ); + + //console.log (relax.initial_grid[0]); + if (relax.has_unclassified) { + relax.initial.unc_mean = ((selection.io.extract_global_MLE_re (relax.final_partitioned_mg_results, "^" + terms.parameters.omega_ratio + ".+`relax.unclassified_branches_name`.+"))["0"])[terms.fit.MLE]; + + relax.initial_grid = utility.Map (relax.initial_grid, "_v_", + 'relax._renormalize_with_weights (_v_, "relax.distribution_uc", relax.initial.unc_mean)' + ); + } + //console.log (relax.initial_grid[0]); + +} //------------------------------------ @@ -997,6 +1013,7 @@ function relax.FitMainTestPair (prompt) { if (relax.do_lhc) { + relax.nm.precision = -0.00025*relax.final_partitioned_mg_results[terms.fit.log_likelihood]; //parameters.DeclareGlobalWithRanges ("relax.bl.scaler", 1, 0, 1000); @@ -1126,11 +1143,12 @@ function relax.FitMainTestPair (prompt) { relax.inferred_distribution_ref2 = parameters.GetStickBreakingDistribution (models.codon.BS_REL.ExtractMixtureDistribution (relax.model_object_map ["relax.reference"])) % 0; selection.io.report_dnds (relax.inferred_distribution_ref2); - + selection.io.json_store_setting (relax.json, "convergence-unstable-alernative", { + {estimators.GetGlobalMLE (relax.alternative_model.fit.take2,terms.relax.k), relax.alternative_model.fit.take2 [terms.fit.log_likelihood]} + }); - if (relax.alternative_model.fit.take2 [terms.fit.log_likelihood] > relax.alternative_model.fit[terms.fit.log_likelihood]) { + if (relax.alternative_model.fit.take2 [terms.fit.log_likelihood] > relax.alternative_model.fit.take2[terms.fit.log_likelihood]) { relax.stash_fitted_K = relax.fitted.K; - io.ReportProgressMessageMD("RELAX", "alt-2", "\n### Potential for highly unreliable K inference due to multiple local maxima in the likelihood function, treat results with caution "); io.ReportProgressMessageMD("RELAX", "alt-2", "> Relaxation parameter reset to opposite mode of evolution from that obtained in the initial optimization."); @@ -1144,7 +1162,6 @@ function relax.FitMainTestPair (prompt) { selection.io.json_store_setting (relax.json, "convergence-unstable", { {relax.fitted.K, relax.alternative_model.fit.take2 [terms.fit.log_likelihood]} {relax.stash_fitted_K, relax.alternative_model.fit [terms.fit.log_likelihood]} - }); io.ReportProgressMessageMD("RELAX", "alt-2", "* The following rate distribution was inferred for **reference** branches"); @@ -1879,12 +1896,14 @@ 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); + mean = Max (mean, 1e-3); over_one = m[d-1][0] * m[d-1][1]; @@ -1907,16 +1926,16 @@ lfunction relax._renormalize_with_weights (v, distro, mean) { under_one = +(m[-1][0] $ m[-1][1]); - for (i = 0; i < d; i+=1) { - m[i][0] = m[i][0] * mean / under_one; + if (under_one > 0) { + 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); + wts = parameters.SetStickBreakingDistributionWeigths (m[-1][1]); + for (i = 0; i < d; i+=1) { (v[((^distro)["rates"])[i]])[^"terms.fit.MLE"] = m[i][0]; diff --git a/res/TemplateBatchFiles/SelectionAnalyses/modules/shared-load-file.bf b/res/TemplateBatchFiles/SelectionAnalyses/modules/shared-load-file.bf index 626fbd277..dcd25bf74 100644 --- a/res/TemplateBatchFiles/SelectionAnalyses/modules/shared-load-file.bf +++ b/res/TemplateBatchFiles/SelectionAnalyses/modules/shared-load-file.bf @@ -56,8 +56,7 @@ function load_nuc_file (prefix) { utility.SetEnvVariable(utility.getGlobalValue ("terms.trees.data_for_neighbor_joining"), None); - - utility.SetEnvVariable(utility.getGlobalValue ("terms.trees.data_for_neighbor_joining"), None); + utility.SetEnvVariable(utility.getGlobalValue ("terms.trees.data_for_neighbor_joining"), None); /** this will return a dictionary of partition strings and trees; one set per partition, as in { @@ -176,8 +175,9 @@ function load_file (prefix) { codon_data_info [utility.getGlobalValue("terms.data.sites")] = {1, file_count}; codon_data_info [utility.getGlobalValue("terms.data.sample_size")] = 0; datasets = {}; - - for (i, filepath; in; file_list) { + + for (i = 0; i < file_count; i+=1) { + filepath = file_list[i]; datasets[i] = prefix+".codon_data_" + i; if (+i == 0) { codon_data_info [filepath] = @@ -472,14 +472,18 @@ function doGTR (prefix) { if (run_gtr) { + utility.ToggleEnvVariable ("TREE_DO_NOT_AUTO_RENAME", TRUE); + gtr_results = estimators.FitGTR(filter_names, trees, gtr_results); - if (Type (save_intermediate_fits) == "AssociativeList") { + utility.ToggleEnvVariable ("TREE_DO_NOT_AUTO_RENAME", None); + + if (Type (save_intermediate_fits) == "AssociativeList") { (save_intermediate_fits[^"terms.data.value"])["GTR"] = gtr_results; io.SpoolJSON (save_intermediate_fits[^"terms.data.value"],save_intermediate_fits[^"terms.data.file"]); - } + } } KeywordArgument ("kill-zero-lengths", "Automatically delete internal zero-length branches for computational efficiency (will not affect results otherwise)", "Yes"); @@ -698,8 +702,12 @@ function doPartitionedMGModel (prefix, keep_lf, model, constraint) { etc */ scaler_variables = utility.PopulateDict (0, partition_count, "`prefix`.scaler_prefix + '_' + _k_", "_k_"); + + for (_value_; in; scaler_variables) { + parameters.DeclareGlobal(_value_, None); + parameters.SetValue(_value_, 3); + } - utility.ForEach (scaler_variables, "_value_", "parameters.DeclareGlobal(_value_, None);parameters.SetValue(_value_, 3);"); run_mg94 = TRUE; @@ -720,6 +728,7 @@ function doPartitionedMGModel (prefix, keep_lf, model, constraint) { } if (run_mg94) { + partitioned_mg_results = estimators.FitCodonModel (filter_names, trees, model, codon_data_info [utility.getGlobalValue("terms.code")], { utility.getGlobalValue("terms.run_options.model_type"): utility.getGlobalValue("terms.local"), utility.getGlobalValue("terms.run_options.proportional_branch_length_scaler"): scaler_variables, diff --git a/res/TemplateBatchFiles/SeqAlignmentCodonShared.ibf b/res/TemplateBatchFiles/SeqAlignmentCodonShared.ibf index f9f21ffc2..24876b7c8 100644 --- a/res/TemplateBatchFiles/SeqAlignmentCodonShared.ibf +++ b/res/TemplateBatchFiles/SeqAlignmentCodonShared.ibf @@ -162,6 +162,7 @@ function cSM2partialSMs(_cdnScoreMatrix, penalties) m3x1[ thisCodon ][ 3 * d1 + 0 ] = max100 - p3x1; m3x1[ thisCodon ][ 3 * d1 + 1 ] = max010 - p3x1; m3x1[ thisCodon ][ 3 * d1 + 2 ] = max001 - p3x1; + } } diff --git a/res/TemplateBatchFiles/libv3/all-terms.bf b/res/TemplateBatchFiles/libv3/all-terms.bf index 5077b4fa9..8bf49d2b7 100644 --- a/res/TemplateBatchFiles/libv3/all-terms.bf +++ b/res/TemplateBatchFiles/libv3/all-terms.bf @@ -17,6 +17,7 @@ namespace terms{ codons = "codons"; codon = "codon"; sense_codons = "sense"; + orf = "ORF"; nucleotide = "nucleotide"; dinucleotide = "dinucleotide"; binary = "binary"; @@ -154,6 +155,9 @@ namespace terms{ function AddCategory (term, categoryID) { return term + " for category " + categoryID; } + function MeanScaler (term) { + return "Mean scaler variable for " + term; + } /* Terms accompanying category defintions */ diff --git a/res/TemplateBatchFiles/libv3/models/codon/MG_REV.bf b/res/TemplateBatchFiles/libv3/models/codon/MG_REV.bf index 783bfb5f0..ed6a68670 100644 --- a/res/TemplateBatchFiles/libv3/models/codon/MG_REV.bf +++ b/res/TemplateBatchFiles/libv3/models/codon/MG_REV.bf @@ -124,8 +124,6 @@ function models.codon.MG_REV._DefineQ(mg_rev, namespace) { */ function models.codon.MG_REV.set_branch_length(model, value, parameter) { - - if (model[terms.model.type] == terms.global) { // boost numeric branch length values by a factor of 3 @@ -180,8 +178,7 @@ function models.codon.MG_REV.set_branch_length(model, value, parameter) { } else { models.codon.MG_REV.set_branch_length.lp = 0; - - + if (parameters.IsIndependent(models.codon.MG_REV.set_branch_length.alpha.p)) { Eval(models.codon.MG_REV.set_branch_length.alpha.p + ":=(" + value[terms.model.branch_length_scaler] + ")*" + value[terms.branch_length]); models.codon.MG_REV.set_branch_length.lp += 1; @@ -194,9 +191,9 @@ function models.codon.MG_REV.set_branch_length(model, value, parameter) { return models.codon.MG_REV.set_branch_length.lp; } } else { - + if (parameters.IsIndependent(models.codon.MG_REV.set_branch_length.alpha.p)) { // alpha is unconstrained; - if (parameters.IsIndependent(models.codon.MG_REV.set_branch_length.beta.p)) { // beta is unconstrained; + if (parameters.IsIndependent(models.codon.MG_REV.set_branch_length.beta.p)) { // beta is unconstrained; models.codon.MG_REV.set_branch_length.lp = parameters.NormalizeRatio(Eval(models.codon.MG_REV.set_branch_length.beta.p), Eval(models.codon.MG_REV.set_branch_length.alpha.p)); parameters.SetConstraint(models.codon.MG_REV.set_branch_length.beta, models.codon.MG_REV.set_branch_length.alpha + "*" + models.codon.MG_REV.set_branch_length.lp, ""); @@ -214,10 +211,11 @@ function models.codon.MG_REV.set_branch_length(model, value, parameter) { parameters.SetConstraint ( models.codon.MG_REV.set_branch_length.alpha.p, models.codon.MG_REV.set_branch_length.alpha, ""); parameters.SetConstraint ( models.codon.MG_REV.set_branch_length.beta, models.codon.MG_REV.set_branch_length.beta.p, ""); - + ConvertBranchLength (models.codon.MG_REV.set_branch_length.lp, model[terms.model.branch_length_string], ^models.codon.MG_REV.set_branch_length.alpha, 3*value); - ExecuteCommands("FindRoot (models.codon.MG_REV.set_branch_length.lp,(" + model[terms.model.branch_length_string] + ")-(" + 3*value + ")," + models.codon.MG_REV.set_branch_length.alpha + ",0,10000)"); + + //ExecuteCommands("FindRoot (models.codon.MG_REV.set_branch_length.lp,(" + model[terms.model.branch_length_string] + ")-(" + 3*value + ")," + models.codon.MG_REV.set_branch_length.alpha + ",0,10000)"); Eval("`models.codon.MG_REV.set_branch_length.alpha.p` =" + models.codon.MG_REV.set_branch_length.lp); diff --git a/res/TemplateBatchFiles/libv3/models/codon/MSS.bf b/res/TemplateBatchFiles/libv3/models/codon/MSS.bf index b77236695..f4100afe8 100644 --- a/res/TemplateBatchFiles/libv3/models/codon/MSS.bf +++ b/res/TemplateBatchFiles/libv3/models/codon/MSS.bf @@ -27,11 +27,10 @@ lfunction model.codon.MSS.prompt_and_define_freq (type, code, freq) { partitioning_option = io.SelectAnOption ( { {"Full", "Each set of codons mapping to the same amino-acid class have a separate substitution rate (Valine == neutral)"} - {"SynREV", "Each set of codons mapping to the same amino-acid class have a separate substitution rate (Valine == neutral)"} - {"SynREVFull", "Each set of codons mapping to the same amino-acid class have a separate substitution rate (no reference)"} + {"SynREV", "Each set of codons mapping to the same amino-acid class have a separate substitution rate (mean = 1)"} {"SynREV2", "Each pair of synonymous codons mapping to the same amino-acid class and separated by a transition have a separate substitution rate (no rate scaling))"} {"SynREV2g", "Each pair of synonymous codons mapping to the same amino-acid class and separated by a transition have a separate substitution rate (Valine == neutral). All between-class synonymous substitutions share a rate."} - {"SynREVCodon", "Each codon pair that is exchangeable gets its own substitution rate (fully estimated)"} + {"SynREVCodon", "Each codon pair that is exchangeable gets its own substitution rate (fully estimated, mean = 1)"} {"Random", "Random partition (specify how many classes; largest class = neutral)"} {"File", "Load a TSV partition from file (prompted for neutral class)"} }, @@ -81,7 +80,11 @@ lfunction model.codon.MSS.prompt_and_define_freq (type, code, freq) { } if (partitioning_option == "SynREV") { return models.codon.MSS.ModelDescription(type, code, - {^"terms.model.MSS.codon_classes" : mapping, ^"terms.model.MSS.neutral" : "V"} + { + ^"terms.model.MSS.codon_classes" : mapping, + ^"terms.model.MSS.normalize" : TRUE + //^"terms.model.MSS.neutral" : "V" + } ); } if (partitioning_option == "SynREVFull") { @@ -90,7 +93,7 @@ lfunction model.codon.MSS.prompt_and_define_freq (type, code, freq) { ); } return models.codon.MSS.ModelDescription(type, code, - {^"terms.model.MSS.codon_classes" : mapping_codon, ^"terms.model.MSS.neutral" : "", ^"terms.model.MSS.normalize" : FALSE} + {^"terms.model.MSS.codon_classes" : mapping_codon, ^"terms.model.MSS.normalize" : TRUE} ); } @@ -210,6 +213,7 @@ lfunction models.codon.MSS.ModelDescription(type, code, codon_classes) { m[^"terms.model.MSS.between"] = codon_classes [^"terms.model.MSS.between"]; } + if (codon_classes[utility.getGlobalValue("terms.model.MSS.normalize")]) { m[utility.getGlobalValue("terms.model.post_definition")] = "models.codon.MSS.post_definition"; } @@ -221,13 +225,24 @@ lfunction models.codon.MSS.ModelDescription(type, code, codon_classes) { lfunction models.codon.MSS.post_definition (model) { // TBD - vars = {}; - weights = {}; - for (i; in; model.GetParameters_RegExp (model,"^" + utility.getGlobalValue ("terms.parameters.synonymous_rate"))) { - vars + i; - weights + 1; + rates = model.GetParameters_RegExp (model,"^" + utility.getGlobalValue ("terms.parameters.synonymous_rate")); + D = utility.Array1D (rates); + w = 1 / D; + + vars = {1,D}; + weights = {1,D}; + terms = {1,D}; + + + i = 0; + for (t,n; in; rates) { + terms[i] = t; + vars[i] = n; + weights[i] = w; + i += 1; } - parameters.ConstrainMeanOfSet (vars, weights, 1, "beavis"); + ((model[utility.getGlobalValue("terms.parameters")])[utility.getGlobalValue("terms.global")]) * (parameters.ConstrainMeanOfSetWithTerms (vars, weights, terms, 1, model[^"terms.id"])[utility.getGlobalValue("terms.global")]); + model = models.generic.post.definition (model); } //---------------------------------------------------------------------------------------------------------------- diff --git a/res/TemplateBatchFiles/libv3/models/model_functions.bf b/res/TemplateBatchFiles/libv3/models/model_functions.bf index 2af115b59..8ec78ab92 100644 --- a/res/TemplateBatchFiles/libv3/models/model_functions.bf +++ b/res/TemplateBatchFiles/libv3/models/model_functions.bf @@ -345,7 +345,7 @@ function model.generic.DefineMixtureModel (model_spec, id, arguments, data_filte } /** - * @name model.generic.GetLocalParameter + * @name models.generic.post.definition * @param {Model} model * @returns {String} */ @@ -400,7 +400,7 @@ function models.generic.ConstrainBranchLength (model, value, parameter) { */ function models.generic.SetBranchLength (model, value, parameter) { - if (Abs((model[terms.parameters])[terms.local]) >= 1) { + if (Abs((model[terms.parameters])[terms.local]) >= 1) { if (Type (model [terms.model.branch_length_string]) == "String") { models.generic.SetBranchLength.expression = model [terms.model.branch_length_string]; @@ -446,11 +446,17 @@ function models.generic.SetBranchLength (model, value, parameter) { } - if (parameters.IsIndependent (models.generic.SetBranchLength.bl.p) || model[terms.model.branch_length_override]) { + if (parameters.IsIndependent (models.generic.SetBranchLength.bl.p) || model[terms.model.branch_length_override]) { if (Type (value) == "AssociativeList") { if (Abs (models.generic.SetBranchLength.expression)) { + //console.log (models.generic.SetBranchLength.expression); + //console.log (value[terms.branch_length]); + //console.log (busted._shared_srv.rv_gdd); + //console.log (^(value[terms.model.branch_length_scaler])); ConvertBranchLength (models.generic.SetBranchLength.t, models.generic.SetBranchLength.expression, ^models.generic.SetBranchLength.bl, value[terms.branch_length]); + //console.log ("**" + Eval (models.generic.SetBranchLength.expression)); + 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); @@ -599,7 +605,7 @@ lfunction models.BindGlobalParametersDeferred (models, filter) { candidate_set[_key_] = 1; } } - + candidate_set = utility.Keys (candidate_set); constraints_set = {}; diff --git a/res/TemplateBatchFiles/libv3/models/parameters.bf b/res/TemplateBatchFiles/libv3/models/parameters.bf index a5c4a6bec..3556f448d 100644 --- a/res/TemplateBatchFiles/libv3/models/parameters.bf +++ b/res/TemplateBatchFiles/libv3/models/parameters.bf @@ -218,6 +218,40 @@ function parameters.SetValues(set) { } } +/** + * Ensures that the mean of parameters in a set is maintained + * @name parameters.ConstrainMeanOfSet + * @param {Dict/Matrix} set - list of variable ids + * @param {Dict/Matrix} weights - weights to apply + * @param {Number} mean - desired mean + * @param {String} namespace - desired mean + * @returns nothing + */ +lfunction parameters.ConstrainMeanOfSetWithTerms (set, weights, terms, mean, namespace) { + + + if (Type (set) == "Matrix") { + unscaled = utility.Map (set, "_name_", "_name_ + '_scaler_variable'"); + constraint = utility.MapWithKey (unscaled, "_key_", "_name_", "_name_ + '*' + `&weights`[_key_[0]+_key_[1]]"); + } + else { + return; + } + + scaler_variables = {}; + utility.ForEach (unscaled, "_name_", 'parameters.DeclareGlobal (_name_, null)'); + global_scaler = namespace + ".scaler_variable"; + parameters.SetConstraint (global_scaler, Join ("+", constraint), "global"); + + for (i,_name_; in; set) { + scaler_variables[terms.MeanScaler(terms[i])] = _name_ + "_scaler_variable"; + parameters.SetValue (_name_ + "_scaler_variable", _name_); + parameters.SetConstraint (_name_, "(" + mean + ")*" + _name_ + "_scaler_variable/`global_scaler`", ""); + } + + return {^'terms.global' : scaler_variables}; +} + /** * Ensures that the mean of parameters in a set is maintained * @name parameters.ConstrainMeanOfSet @@ -241,16 +275,18 @@ lfunction parameters.ConstrainMeanOfSet (set, weights, mean, namespace) { } } + scaler_variables = {}; utility.ForEach (unscaled, "_name_", 'parameters.DeclareGlobal (_name_, null)'); global_scaler = namespace + ".scaler_variable"; parameters.SetConstraint (global_scaler, Join ("+", constraint), "global"); - utility.ForEach (set, "_name_", ' - `&scaler_variables`["Mean scaler variable for " + _name_] = _name_ + "_scaler_variable"; - parameters.SetValue (_name_ + "_scaler_variable", _name_); - parameters.SetConstraint (_name_, "(" + `&mean` + ")*" + _name_ + "_scaler_variable/`global_scaler`", ""); - '); + for (_name_; in; set) { + scaler_variables[terms.MeanScaler(_name_)] = _name_ + "_scaler_variable"; + parameters.SetValue (_name_ + "_scaler_variable", _name_); + parameters.SetConstraint (_name_, "(" + mean + ")*" + _name_ + "_scaler_variable/`global_scaler`", ""); + } + return {^'terms.global' : scaler_variables}; } @@ -743,11 +779,14 @@ lfunction parameters.SetStickBreakingDistributionWeigths (weights) { left_over = (1-w[0]); for (i = 1; i < rate_count - 1; i += 1) { - w [i] = weights[i] / left_over; + if (left_over > 0) { + w [i] = weights[i] / left_over; + } else { + w [i] = 0; + } left_over = left_over * (1-w[i]); } - //w[i] = left_over; return w; } diff --git a/res/TemplateBatchFiles/libv3/tasks/alignments.bf b/res/TemplateBatchFiles/libv3/tasks/alignments.bf index 0d4c0d3a9..bd477b7b5 100644 --- a/res/TemplateBatchFiles/libv3/tasks/alignments.bf +++ b/res/TemplateBatchFiles/libv3/tasks/alignments.bf @@ -732,8 +732,18 @@ lfunction alignments.TranslateCodonsToAminoAcidsWithAmbiguities (sequence, offse }; upper_bound = Abs (try_run); + orf = {{0,-1}}; + last_nonstop = 0; + last_stop = -1; for (i = 0; i < upper_bound; i+=1) { if (try_run[i] / "X") { // has_stop + if (i > last_nonstop) { + if (i-last_stop-1 > orf[1]-orf[0]) { + orf[0] = last_stop+1; + orf[1] = i-1; + } + } + last_stop = i; translation * "X"; frame_result [^"terms.stop_codons"] += 1; if (i == upper_bound - 1) { @@ -742,9 +752,16 @@ lfunction alignments.TranslateCodonsToAminoAcidsWithAmbiguities (sequence, offse } else { translation * (try_run[i])["INDEXORDER"][0]; frame_result [^"terms.sense_codons"] += 1; + last_nonstop = i; } } - + + if (last_nonstop-last_stop-1 > orf[1]-orf[0]) { + orf[0] = last_stop+1; + orf[1] = last_nonstop; + } + + frame_result [^"terms.orf"] = orf; translation * 0; frame_result [utility.getGlobalValue ("terms.data.sequence")] = translation; diff --git a/res/TemplateBatchFiles/libv3/tasks/estimators.bf b/res/TemplateBatchFiles/libv3/tasks/estimators.bf index 992ed68a7..ba26ebf4b 100644 --- a/res/TemplateBatchFiles/libv3/tasks/estimators.bf +++ b/res/TemplateBatchFiles/libv3/tasks/estimators.bf @@ -547,6 +547,30 @@ function estimators.ApplyExistingEstimatesToTree (_tree_name, model_descriptions return estimators.ApplyExistingEstimatesToTree.constraint_count; } +/** + + Apply existing constraints and override previously constrained BL estimates + + * @name estimators.ApplyExistingEstimatesOverride + * @param {String} likelihood_function_id + * @param {Dictionary} model_descriptions + * @param {Matrix} initial_values + * @param branch_length_conditions + * @returns estimators.ApplyExistingEstimates.df_correction - Abs(estimators.ApplyExistingEstimates.keep_track_of_proportional_scalers); + */ + + +function estimators.ApplyExistingEstimatesOverride(likelihood_function_id, model_descriptions, initial_values, branch_length_conditions) { + for (i, v; in; model_descriptions) { + v[^"terms.model.branch_length_override"] = TRUE; + } + estimators.ApplyExistingEstimates(likelihood_function_id, model_descriptions, initial_values, branch_length_conditions); + for (i, v; in; model_descriptions) { + v - ^"terms.model.branch_length_override"; + } +} + + /** * @name * @param {String} likelihood_function_id @@ -657,7 +681,7 @@ lfunction estimators.FitExistingLF (lf_id, model_objects) { /** * Makes a likelihood function object with the desired parameters - * @name estimators.FitLF + * @name estimators.BuildLFObject * @param {Matrix} data_filters_list - a vector of {DataFilter}s * @param {Matrix} tree_list - a vector of {Tree}s * @param model_map @@ -776,17 +800,20 @@ lfunction estimators.FitLF(data_filter, tree, model_map, initial_values, model_o can_do_restarts = null; if (utility.Has (run_options, utility.getGlobalValue("terms.search_grid"),"AssociativeList")) { + 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"); + //console.log (run_options [utility.getGlobalValue("terms.search_grid")]); + //console.log (grid_results); + 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]]; @@ -796,11 +823,9 @@ lfunction estimators.FitLF(data_filter, tree, model_map, initial_values, model_o if (null == can_do_restarts) { best_value = Max (grid_results, 1); parameters.SetValues ((run_options [utility.getGlobalValue("terms.search_grid")])[best_value["key"]]); + //console.log ((run_options [utility.getGlobalValue("terms.search_grid")])[best_value["key"]]); + estimators.ApplyExistingEstimatesOverride (&likelihoodFunction, model_objects, initial_values, run_options[utility.getGlobalValue("terms.run_options.proportional_branch_length_scaler")]); } - //console.log (best_value); - //console.log ((run_options [utility.getGlobalValue("terms.search_grid")])[best_value["key"]]); - //console.log (can_do_restarts); - //assert (0); } @@ -810,14 +835,18 @@ lfunction estimators.FitLF(data_filter, tree, model_map, initial_values, model_o utility.ToggleEnvVariable("PRODUCE_OPTIMIZATION_LOG", 1); } + //utility.ToggleEnvVariable("VERBOSITY_LEVEL", 101); + //console.log (run_options[utility.getGlobalValue("terms.run_options.optimization_settings")]); if (Type (can_do_restarts) == "AssociativeList") { io.ReportProgressBar("", "Working on crude initial optimizations"); bestlog = -1e100; for (i = 0; i < Abs (can_do_restarts); i += 1) { - parameters.SetValues (can_do_restarts[i]); + + parameters.SetValues (can_do_restarts[i]); + estimators.ApplyExistingEstimatesOverride (&likelihoodFunction, model_objects, initial_values, run_options[utility.getGlobalValue("terms.run_options.proportional_branch_length_scaler")]); if (utility.Has (run_options,utility.getGlobalValue("terms.run_options.optimization_settings"),"AssociativeList")) { - Optimize (mles, likelihoodFunction, run_options[utility.getGlobalValue("terms.run_options.optimization_settings")]); + Optimize (mles, likelihoodFunction, run_options[utility.getGlobalValue("terms.run_options.optimization_settings")]); } else { Optimize (mles, likelihoodFunction); } @@ -841,6 +870,8 @@ lfunction estimators.FitLF(data_filter, tree, model_map, initial_values, model_o results = estimators.ExtractMLEs( & likelihoodFunction, model_objects); results[utility.getGlobalValue ("terms.fit.log_likelihood")] = mles[1][0]; } + + if (optimization_log) { @@ -1178,7 +1209,11 @@ lfunction estimators.FitCodonModel(codon_data, tree, generator, genetic_code, op Assumes that option["partitioned-omega"] is a dictionary where each partition has an entry (0-index based), which itself is a dictionary of the form: "branch-name" : "branch-set" */ - utility.ForEach(option[utility.getGlobalValue("terms.run_options.partitioned_omega")], "_value_", "utility.AddToSet(`&partition_omega`,utility.UniqueValues(_value_))"); + + + for (_value_; in; option[utility.getGlobalValue("terms.run_options.partitioned_omega")]) { + utility.AddToSet(partition_omega,utility.UniqueValues(_value_)); + } } @@ -1192,10 +1227,16 @@ lfunction estimators.FitCodonModel(codon_data, tree, generator, genetic_code, op new_globals = {}; - utility.ForEachPair(partition_omega, "_key_", "_value_", - '`&new_globals` [_key_] = (`&name_space` + ".omega_" + Abs (`&new_globals`)); model.generic.AddGlobal (`&mg_rev`, `&new_globals` [_key_] , (utility.getGlobalValue("terms.parameters.omega_ratio")) + " for *" + _key_ + "*")'); + + for (_key_, _value_; in; partition_omega) { + new_globals[_key_] = name_space + ".omega_" + Abs (new_globals); + model.generic.AddGlobal (mg_rev, + new_globals[_key_], + (utility.getGlobalValue("terms.parameters.omega_ratio")) + " for *" + _key_ + "*" + ); + } parameters.DeclareGlobal(new_globals, None); - + /** now replicate the local constraint for individual branches @@ -1204,12 +1245,13 @@ lfunction estimators.FitCodonModel(codon_data, tree, generator, genetic_code, op alpha = model.generic.GetLocalParameter(mg_rev, utility.getGlobalValue("terms.parameters.synonymous_rate")); beta = model.generic.GetLocalParameter(mg_rev, utility.getGlobalValue("terms.parameters.nonsynonymous_rate")); + io.CheckAssertion("None!=`&alpha` && None!=`&beta`", "Could not find expected local synonymous and non-synonymous rate parameters in \`estimators.FitMGREV\`"); SetParameter (DEFER_CONSTRAINT_APPLICATION, 1, 0); apply_constraint: = component_tree + "." + node_name + "." + beta + ":=" + component_tree + "." + node_name + "." + alpha + "*" + new_globals[branch_map[node_name]]; - + for (i = 0; i < components; i += 1) { component_tree = lf_components[2 * i + 1]; ClearConstraints( * component_tree); diff --git a/res/TemplateBatchFiles/libv3/tasks/mpi.bf b/res/TemplateBatchFiles/libv3/tasks/mpi.bf index a94a61d01..b544c20a7 100644 --- a/res/TemplateBatchFiles/libv3/tasks/mpi.bf +++ b/res/TemplateBatchFiles/libv3/tasks/mpi.bf @@ -347,8 +347,12 @@ namespace mpi { for (i = 0; i < task_count; i+=1) { parameters.SetValues (tasks[task_ids[i]]); + //console.log (tasks[task_ids[i]]); + //console.log (values[1]); estimators.ApplyExistingEstimates (lf_id, values[0], values[1], values[2]); + //GetInformation (i, ^"busted._shared_srv.rv_gdd"); //fprintf (stdout, ^lf_id); + //assert (0); LFCompute (^lf_id, ll); results [task_ids[i]] = ll; io.ReportProgressBar("", "Grid result " + i + " = " + ll + " (best = " + Max (results, 0)[^"terms.data.value"] + ")"); diff --git a/res/TemplateBatchFiles/libv3/tasks/trees.bf b/res/TemplateBatchFiles/libv3/tasks/trees.bf index 9a2b09fdf..33289b3ab 100644 --- a/res/TemplateBatchFiles/libv3/tasks/trees.bf +++ b/res/TemplateBatchFiles/libv3/tasks/trees.bf @@ -312,6 +312,7 @@ lfunction trees.LoadAnnotatedTreeTopology.match_partitions(partitions, mapping) partition_count = Rows(partitions); partrees = {}; + tree_matrix = utility.GetEnvVariable("NEXUS_FILE_TREE_MATRIX"); diff --git a/src/core/alignment.cpp b/src/core/alignment.cpp index bda824642..325a6ad38 100644 --- a/src/core/alignment.cpp +++ b/src/core/alignment.cpp @@ -563,13 +563,14 @@ inline void BacktrackAlignCodon( signed char * const edit_ops 2. there's also an opportinity to use bitmasks */ + switch ( code ) { // match case HY_111_111: r -= 3; q -= 3; - edit_ops [edit_ptr] = 1; - edit_ops [edit_ptr+1] = 1; - edit_ops [edit_ptr+2] = 1; + edit_ops [edit_ptr] = 0; + edit_ops [edit_ptr+1] = 0; + edit_ops [edit_ptr+2] = 0; edit_ptr += 3; return; @@ -594,14 +595,13 @@ inline void BacktrackAlignCodon( signed char * const edit_ops unsigned char r_str[ 5 ] = { 1, 1, 1, 1, 1 }, // which characters are we taking from the reference (up to 5); 1 take, 0 leave - q_str[ 5 ] = { 1, 1, 1, 1, 1 }, + q_str[ 5 ] = { 1, 1, 1, 1, 1 }; // which characters are we taking from the query (up to 5) - idx = 2; + int idx = 2; // how many characters are we taking (max over both strings); index so the actual number is +1 // can be 2,3, or 4 - switch (code) { // 3x2 case HY_111_110: @@ -1146,7 +1146,7 @@ double AlignStrings( char const * r_str // fill in the edit_ops with the difference // between q_len and j - for (long k = index_R; k < q_len; ++k ) { + for (long k = index_Q; k < q_len; ++k ) { edit_ops[ edit_ptr++ ] = 1; } } @@ -1336,10 +1336,13 @@ double AlignStrings( char const * r_str } + if (!took_local_shortcut) { // for anything that remains, // don't forget it!!! // reference + + while ( --index_R >= 0 ) edit_ops[ edit_ptr++ ] = -1; @@ -1764,6 +1767,7 @@ double LinearSpaceAlign ( const char *s1, // first string char *ha ) { + if (to2 == from2 || to1 == from1) { return 0; } diff --git a/src/core/batchlanruntime.cpp b/src/core/batchlanruntime.cpp index 50e5591ea..c14acbb68 100644 --- a/src/core/batchlanruntime.cpp +++ b/src/core/batchlanruntime.cpp @@ -1158,7 +1158,7 @@ bool _ElementaryCommand::HandleAlignSequences(_ExecutionList& current_progr _SimpleList ops (reference_sequence->length()+2UL,-2,0); ops.list_data[reference_sequence->length()+1] = sequence2->length(); ops.list_data[0] = -1; - + score = LinearSpaceAlign (reference_sequence->get_str(), sequence2->get_str(), reference_sequence->length(), @@ -1181,7 +1181,7 @@ bool _ElementaryCommand::HandleAlignSequences(_ExecutionList& current_progr data_buffers, 0, alignment_route); - + delete[] alignment_route; _StringBuffer *result1 = new _StringBuffer (reference_sequence->length() + 1UL), diff --git a/src/core/formula.cpp b/src/core/formula.cpp index 721eb74a0..88c7daa90 100644 --- a/src/core/formula.cpp +++ b/src/core/formula.cpp @@ -1933,6 +1933,12 @@ _Variable * _Formula::Dereference (bool ignore_context, _hyExecutionContext* the HBLObjectRef _Formula::Compute (long startAt, _VariableContainer const * nameSpace, _List* additionalCacheArguments, _String* errMsg, long valid_type, bool can_cache, bool skip_comparison) { // compute the value of the formula // TODO SLKP 20170925 Needs code review + + if (simpleExpressionStatus) { + HandleApplicationError("Internal error, calling _Formula::Compute on a compiled _Formula object"); + } + + _Stack * scrap_here; current_formula_being_computed = this; if (theFormula.empty()) { @@ -2367,16 +2373,17 @@ bool _Formula::ConvertToSimple (_AVLList& variable_index) { //__________________________________________________________________________________ void _Formula::ConvertFromSimple (_AVLList const& variableIndex) { - delete [] simpleExpressionStatus; - simpleExpressionStatus = nil; ConvertFromSimpleList (*variableIndex.dataList); } //__________________________________________________________________________________ void _Formula::ConvertFromSimpleList (_SimpleList const& variableIndex) { - if (theFormula.empty()) { - return; - } + delete [] simpleExpressionStatus; + simpleExpressionStatus = nil; + + if (theFormula.empty()) { + return; + } for (unsigned long i=0UL; itheNumber) { @@ -2923,7 +2944,7 @@ long _Formula::ObjectClass (void) { if (theStack.StackDepth()) { return ((HBLObjectRef)theStack.theStack.list_data[0])->ObjectClass(); } - + HBLObjectRef res = Compute(); if (res) { @@ -3034,3 +3055,27 @@ void _Formula::ConvertFromTree (void) { } } +//__________________________________________________________________________________ +bool _Formula::ProcessFormulaForConverstionToSimple (long& stack_length, + _AVLList& var_list, + _SimpleList& new_formulas, + _SimpleList& references, + _AVLListX& formula_strings, + bool run_all) { + + if (run_all || AmISimple(stack_length,var_list)) { + _String * formula_string = (_String*)toStr(kFormulaStringConversionNormal, nil,true); + long fref = formula_strings.Insert(formula_string,new_formulas.lLength); + if (fref < 0) { + references << formula_strings.GetXtra (-fref-1); + DeleteObject (formula_string); + } else { + new_formulas << (long)this; + references << fref; + } + + } else { + return false; + } + return true; +} diff --git a/src/core/fstring.cpp b/src/core/fstring.cpp index 90553c903..3e54070a7 100644 --- a/src/core/fstring.cpp +++ b/src/core/fstring.cpp @@ -598,7 +598,8 @@ HBLObjectRef _FString::SubstituteAndSimplify(HBLObjectRef arguments, HBLObjectRe evaluator.SimplifyConstants(); } - return _returnStringOrUseCache(((_String*)evaluator.toStr(kFormulaStringConversionNormal)), cache); + _String * simplified = ((_String*)evaluator.toStr(kFormulaStringConversionNormal)); + return _returnStringOrUseCache(simplified, cache); } } return new _MathObject; diff --git a/src/core/global_things.cpp b/src/core/global_things.cpp index 58d6d5911..7654e7b02 100644 --- a/src/core/global_things.cpp +++ b/src/core/global_things.cpp @@ -122,7 +122,7 @@ namespace hy_global { kErrorStringMatrixExportError ("Export matrix called with a non-polynomial matrix argument"), kErrorStringNullOperand ("Attempting to operate on an undefined value; this is probably the result of an earlier 'soft' error condition"), kErrorNumerical ("To treat numerical errors as warnings, please specify \"ENV=TOLERATE_NUMERICAL_ERRORS=1;\" as the command line argument. This often resolves the issue, which is indicative of numerical instability."), - kHyPhyVersion = _String ("2.5.60"), + kHyPhyVersion = _String ("2.5.61"), kNoneToken = "None", kNullToken = "null", diff --git a/src/core/include/defines.h b/src/core/include/defines.h index 18fbe3225..9859cc1db 100644 --- a/src/core/include/defines.h +++ b/src/core/include/defines.h @@ -103,10 +103,12 @@ SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. // NeedToExponentiate #define HY_VARIABLE_NOTSET 0x0080 #define HY_VARIABLE_COMPUTING 0x0100 +#define HY_VARIABLE_SET_EXTERNALLY 0x0200 #define HY_VARIABLE_SET 0xFF7F #define HY_DEP_V_INSPECTED_CLR 0xFFF7 #define HY_VARIABLE_COMPUTING_CLR 0xFEFF +#define HY_VARIABLE_SET_EXTERNALLY_CLR 0xFDFF #define HY_VC_CLR_NO_CHECK 0xFFBF #define HY_DEP_CLEAR_MASK 0xFFC7 #define HY_HY_VARIABLE_CHANGED_CLEAR 0xFFFC diff --git a/src/core/include/formula.h b/src/core/include/formula.h index 56c7821d8..f1b03e1dc 100644 --- a/src/core/include/formula.h +++ b/src/core/include/formula.h @@ -205,7 +205,7 @@ class _Formula { */ - bool AmISimple (long& stack_depth, _AVLList& variable_index); + bool AmISimple (long& stack_depth, _AVLList& variable_index) ; long StackDepth (long start_at = 0L, long end_at = -1L) const; /** starting at operation 'start_at', counting up to 'end_at' (-1 == the end), @@ -249,6 +249,25 @@ class _Formula { void ScanFormulaForHBLFunctions (_AVLListX& collection , bool recursive, bool simplify = true); + /** + Process a formula as a part of a batch to convert to simple formulas + + @author SLKP + @date 20240305 + @param stackLength : keep track of the deepest stack for the formulas in the batch + @param var_list : the union of all independent variables that this batch depends on + @param new_formulas : a list of references to unique formulas in the batch + @param references : the index of the unique formula to which this formula refers + @param formula_strings : keep track of formula strings to identify unique ones + @param run_all : if true, process all formulas (not just those which return TRUE from AmISimple) + + */ + bool ProcessFormulaForConverstionToSimple (long& stack_length, + _AVLList& var_list, + _SimpleList& new_formulas, + _SimpleList& references, + _AVLListX& formula_strings, + bool run_all); /** A compute and forget utility function. Parse an expression, optionally check to see that it's of the right type and return the value diff --git a/src/core/include/likefunc.h b/src/core/include/likefunc.h index bb5d7f124..2bccf10ab 100644 --- a/src/core/include/likefunc.h +++ b/src/core/include/likefunc.h @@ -538,6 +538,14 @@ class _LikelihoodFunction: public BaseObj { partition */ + void CompileConstraints (void); + void UncompileConstraints (void); + + /** + 20240305 SLKP + These two functions populate / depopulate 'compiled_constraints' + + */ _List* RecoverAncestralSequencesMarginal (long, _Matrix&,_List const&, bool = false); @@ -721,6 +729,12 @@ class _LikelihoodFunction: public BaseObj { ; _AssociativeList *optimizatonHistory; + + _CompiledMatrixData *compiled_constraints; + /** + SLKP: 20240305 + Use this object for storing compiled + */ #ifdef _OPENMP long lfThreadCount; diff --git a/src/core/include/matrix.h b/src/core/include/matrix.h index c2c677bbe..8b4bd5666 100644 --- a/src/core/include/matrix.h +++ b/src/core/include/matrix.h @@ -72,11 +72,11 @@ struct _CompiledMatrixData { _SimpleFormulaDatum * theStack, * varValues; - hyFloat * formulaValues; + hyFloat * formulaValues; - long * formulaRefs; - long stackDepth; - bool has_volatile_entries; + long * formulaRefs; + long stackDepth; + bool has_volatile_entries; _SimpleList varIndex, formulasToEval; diff --git a/src/core/include/topology.h b/src/core/include/topology.h index 27f367dc2..89905eeeb 100644 --- a/src/core/include/topology.h +++ b/src/core/include/topology.h @@ -63,6 +63,7 @@ struct _TreeTopologyParseSettings { auto_convert_lengths = false; accept_user_lengths = true; ingore_user_inode_names = false; + auto_rename_nodes = true; parser_namespace = kEmptyString; parser_cache = nil; } @@ -82,9 +83,11 @@ struct _TreeTopologyParseSettings { _String inode_prefix, parser_namespace; + bool auto_convert_lengths, accept_user_lengths, - ingore_user_inode_names; + ingore_user_inode_names, + auto_rename_nodes; _AVLListX * parser_cache; }; diff --git a/src/core/include/variable.h b/src/core/include/variable.h index 2ca5f24f4..b845d6eba 100644 --- a/src/core/include/variable.h +++ b/src/core/include/variable.h @@ -100,7 +100,7 @@ class _Variable : public _Constant { if (varValue) { DeleteObject (varValue); varValue = nil;} } - const _Formula * get_constraint (void) const { + _Formula * get_constraint (void) { return varFormula; } diff --git a/src/core/likefunc.cpp b/src/core/likefunc.cpp index b94c9b95e..6f960ec1d 100644 --- a/src/core/likefunc.cpp +++ b/src/core/likefunc.cpp @@ -602,6 +602,7 @@ void _LikelihoodFunction::Init (void) { branchCaches = nil; parameterValuesAndRanges = nil; optimizatonHistory = nil; + compiled_constraints = nil; _variables_changed_during_last_compute = nil; variables_changed_during_last_compute = nil; @@ -778,6 +779,7 @@ bool _LikelihoodFunction::MapTreeTipsToData (long f, _String *errorMessage, b } df->SetMap(tipMatches); } + ReportWarning (_String ("The tips of the tree:") & *t->GetName() &" were matched with the species names from the data in the following numeric order (0-based) "& _String ((_String*)tipMatches.toStr())); } else { _String warnMsg = _String ("The leaf of the tree ") & *t->GetName() &" labeled " &(*(_String*)tips(j)).Enquote() @@ -1933,45 +1935,94 @@ bool _LikelihoodFunction::SendOffToMPI (long index) { //_______________________________________________________________________________________ bool _LikelihoodFunction::PreCompute (void) { - - useGlobalUpdateFlag = true; - // mod 20060125 to only update large globals once - unsigned long i; - - _SimpleList * arrayToCheck = nonConstantDep?nonConstantDep:&indexDep; - - for (i = 0UL; i < arrayToCheck->lLength; i++) { - _Variable* cornholio = LocateVar(arrayToCheck->list_data[i]); - hyFloat tp = cornholio->Compute()->Value(); - if (!cornholio->IsValueInBounds(tp)){ - ReportWarning (_String ("Failing bound checks on ") & *cornholio->GetName() & " = " & _String (tp, "%25.16g")); - //break; + if (compiled_constraints) { + // populate all the independent variables + for (long i = 0; i < compiled_constraints->varIndex.lLength; i++) { + compiled_constraints->varValues[i].value = LocateVar(compiled_constraints->varIndex[i])->Value(); + //printf ("\nArgument %s = %g\n", LocateVar(compiled_constraints->varIndex[i])->GetName()->get_str(), compiled_constraints->varValues[i].value); } - } - - useGlobalUpdateFlag = false; - // mod 20060125 to only update large globals once - - for (unsigned long j=0UL; j < i; j++) { - _Variable* cornholio = LocateVar(arrayToCheck->list_data[j]); - if (cornholio->varFlags&HY_DEP_V_COMPUTED) { - cornholio->varFlags -= HY_DEP_V_COMPUTED; + + // compute all the unique expressions + for (long i = 0; i < compiled_constraints->formulasToEval.lLength; i++) { + compiled_constraints->formulaValues[i] = ((_Formula*)(compiled_constraints->formulasToEval.get (i)))->ComputeSimple(compiled_constraints->theStack, compiled_constraints->varValues); + + //printf ("\nComputed formula %d to %g\n", i, compiled_constraints->formulaValues[i]); + } + + // check the ranges of all the dependent variables + for (long i = 0; i < indexDep.lLength; i++) { + _Variable *ith_dep = GetIthDependentVar(i); + hyFloat dep_value = compiled_constraints->formulaValues[compiled_constraints->formulaRefs[i]]; + if (!ith_dep->IsValueInBounds(dep_value)) { + ReportWarning (_String ("Failing bound checks on ") & *ith_dep->GetName() & " = " & _String (dep_value, "%25.16g")); + } + //printf ("\nSetting %s to %g\n", ith_dep->GetName()->get_str(), dep_value); + ith_dep->SetValue(dep_value); } + return true; + + } else { + useGlobalUpdateFlag = true; + // mod 20060125 to only update large globals once + unsigned long i; + + _SimpleList * arrayToCheck = nonConstantDep?nonConstantDep:&indexDep; + + for (i = 0UL; i < arrayToCheck->lLength; i++) { + _Variable* cornholio = LocateVar(arrayToCheck->list_data[i]); + hyFloat tp = cornholio->Compute()->Value(); + if (!cornholio->IsValueInBounds(tp)){ + ReportWarning (_String ("Failing bound checks on ") & *cornholio->GetName() & " = " & _String (tp, "%25.16g")); + //break; + } + } + + useGlobalUpdateFlag = false; + // mod 20060125 to only update large globals once + + for (unsigned long j=0UL; j < i; j++) { + _Variable* cornholio = LocateVar(arrayToCheck->list_data[j]); + if (cornholio->varFlags&HY_DEP_V_COMPUTED) { + cornholio->varFlags -= HY_DEP_V_COMPUTED; + } + } + + + return (i==arrayToCheck->lLength); } - - - return (i==arrayToCheck->lLength); } //_______________________________________________________________________________________ void _LikelihoodFunction::PostCompute (void) { - _SimpleList * arrayToCheck = nonConstantDep?nonConstantDep:&indexDep; - - //useGlobalUpdateFlag = true; - for (unsigned long i=0; ilLength; i++) { - LocateVar (arrayToCheck->list_data[i])->Compute(); + + if (compiled_constraints) { + // populate all the independent variables + for (long i = 0; i < compiled_constraints->varIndex.lLength; i++) { + compiled_constraints->varValues[i].value = LocateVar(compiled_constraints->varIndex[i])->Value(); + } + + // compute all the unique expressions + for (long i = 0; i < compiled_constraints->formulasToEval.lLength; i++) { + compiled_constraints->formulaValues[i] = ((_Formula*)(compiled_constraints->formulasToEval.get (i)))->ComputeSimple(compiled_constraints->theStack, compiled_constraints->varValues); + } + + // check the ranges of all the dependent variables + for (long i = 0; i < indexDep.lLength; i++) { + _Variable *ith_dep = GetIthDependentVar(i); + hyFloat dep_value = compiled_constraints->formulaValues[compiled_constraints->formulaRefs[i]]; + ith_dep->SetValue(dep_value); + } + + } else { + + _SimpleList * arrayToCheck = nonConstantDep?nonConstantDep:&indexDep; + + //useGlobalUpdateFlag = true; + for (unsigned long i=0; ilLength; i++) { + LocateVar (arrayToCheck->list_data[i])->Compute(); + } } //useGlobalUpdateFlag = false; // mod 20060125 comment out the compute loop; seems redundant @@ -3883,6 +3934,98 @@ void _LikelihoodFunction::LoggerSingleVariable (unsigned long inde LoggerLogL (logL); *((_Vector*) (((_AssociativeList*)this->optimizatonHistory->GetByKey("Parameters")))->GetByKey(*GetIthIndependentName(index))) << GetIthIndependent(index); } +} + +//_______________________________________________________________________________________________ + +void _LikelihoodFunction::UncompileConstraints (void) { + if (!compiled_constraints) return; // nothing to do + + for (unsigned long k = 0; k < compiled_constraints->formulasToEval.lLength; k++) { + //delete ((_Formula*)compiled_constraints->formulasToEval.get(k)); + ((_Formula*)compiled_constraints->formulasToEval.get(k))->ConvertFromSimpleList (compiled_constraints->varIndex); + } + + free (compiled_constraints->theStack); + free (compiled_constraints->varValues); + free (compiled_constraints->formulaRefs); + delete [] compiled_constraints->formulaValues; + + delete compiled_constraints; + compiled_constraints = nil; +} + + +//_______________________________________________________________________________________________ + +void _LikelihoodFunction::CompileConstraints (void) { + if (compiled_constraints) return; // don't REpopulate + + return; + long stackLength = 0L; + + compiled_constraints = new _CompiledMatrixData; + + _SimpleList references ; + + _List flaStringsL; + _AVLListX flaStrings(&flaStringsL); + + _AVLList varList (&compiled_constraints->varIndex); + + bool is_good = true; + for (long i = 0L; i < indexDep.countitems(); i++) { + _Variable * ith_dep = GetIthDependentVar(i); + if (!ith_dep->get_constraint()->ProcessFormulaForConverstionToSimple(stackLength, varList, compiled_constraints->formulasToEval, references, flaStrings, false)) { + is_good = false; + } + } + + // SLKP 20240305 : at this time, if there are any dependent variables in the list of argument variables, bail + + for (long k = 0L; k < varList.countitems(); k++) { + if (!LocateVar(compiled_constraints->varIndex.get (k))->IsIndependent()) { + printf ("\nBAILING -- DEPENDENT 'INDEPENDENT' (%s)\n", LocateVar(compiled_constraints->varIndex.get (k))->theName->get_str()); + is_good = false; + break; + } + } + + + if (!is_good) { + delete compiled_constraints; + compiled_constraints = nil; + return; + } + + compiled_constraints->has_volatile_entries = false; + compiled_constraints->stackDepth = stackLength; + + for (unsigned long k = 0; k < compiled_constraints->formulasToEval.lLength; k++) { + _Formula * ref = (_Formula*)compiled_constraints->formulasToEval.get(k); + compiled_constraints->has_volatile_entries = ref->ConvertToSimple(varList) || compiled_constraints->has_volatile_entries; + } + + compiled_constraints->theStack = (_SimpleFormulaDatum*)MemAllocate (stackLength*sizeof(_SimpleFormulaDatum)); + compiled_constraints->varValues = (_SimpleFormulaDatum*)MemAllocate ((compiled_constraints->varIndex.countitems()>0?varList.countitems():1)*sizeof(_SimpleFormulaDatum)); + long allocation_size = MAX (references.lLength, 1) * sizeof (long); + compiled_constraints->formulaRefs = (long*)MemAllocate (allocation_size); + memcpy (compiled_constraints->formulaRefs, references.list_data, allocation_size); + compiled_constraints->formulaValues = new hyFloat [compiled_constraints->formulasToEval.lLength]; + + /*for (long k = 0L; k < varList.countitems(); k++) { + printf ("Independent %d (%d) => %s\n", k, varListAux.get (k), LocateVar(varListAux.get (k))->theName->get_str()); + }*/ + + printf ( +"\nProcessed LF constraints. \ +\n\tStack depth: %ld\ +\n\tIndependent variables: %ld\ +\n\tUnique expressions: %ld\n", stackLength, varList.countitems(), compiled_constraints->formulasToEval.countitems()); + + + + } //_______________________________________________________________________________________________ @@ -4112,7 +4255,6 @@ _Matrix* _LikelihoodFunction::Optimize (_AssociativeList const * options) SetupLFCaches (); SetupCategoryCaches (); - #ifdef __HYPHYMPI__ if (hy_mpi_node_rank == 0) { @@ -4134,6 +4276,7 @@ _Matrix* _LikelihoodFunction::Optimize (_AssociativeList const * options) } #endif + CompileConstraints (); computationalResults.Clear(); #if !defined __UNIX__ || defined __HEADLESS__ || defined __HYPHYQT__ || defined __HYPHY_GTK__ @@ -4168,6 +4311,7 @@ _Matrix* _LikelihoodFunction::Optimize (_AssociativeList const * options) if (custom_convergence_callback >= 0) { if (GetBFFunctionArgumentCount (custom_convergence_callback) != 2L) { HandleApplicationError("Custom convergence criterion convergence function must have exactly two arguments: current log-L, and an dictionary with id -> value mapping"); + UncompileConstraints (); return new _Matrix; } } @@ -4254,7 +4398,10 @@ _Matrix* _LikelihoodFunction::Optimize (_AssociativeList const * options) precision = get_optimization_setting (kOptimizationPrecision, 0.001); maxItersPerVar = get_optimization_setting (kMaximumIterationsPerVariable, 5000.); - + if (CheckEqual (precision, 0.0)) { + ReportWarning (kOptimizationPrecision & " was set to 0. Resetting to a default value of 0.001"); + precision = 0.001; + } _FString * method_string = get_optimization_setting_string (kOptimizationMethod); @@ -4403,9 +4550,9 @@ _Matrix* _LikelihoodFunction::Optimize (_AssociativeList const * options) hyFloat grad_precision; if (maxSoFar > -INFINITY) { - grad_precision = MIN (1000., -maxSoFar / 100.); + grad_precision = MIN (1000., -maxSoFar / 500.); } else { - grad_precision = -0.01; + grad_precision = -0.005; } if (gradientBlocks.nonempty()) { @@ -4428,7 +4575,7 @@ _Matrix* _LikelihoodFunction::Optimize (_AssociativeList const * options) hyFloat current_precision = MAX(1., precision); while (current_precision > precision) { ConjugateGradientDescent(current_precision, bestSoFar, true); - current_precision *= 0.1; + current_precision *= sqrt(0.1); } ConjugateGradientDescent(precision, bestSoFar, true); } @@ -4452,7 +4599,6 @@ _Matrix* _LikelihoodFunction::Optimize (_AssociativeList const * options) currentPrecision = kOptimizationGradientDescent==7?sqrt(precision):intermediateP; } - if (optimization_mode == kOptimizationCoordinateWise) { bool forward = false; @@ -4730,7 +4876,7 @@ _Matrix* _LikelihoodFunction::Optimize (_AssociativeList const * options) _Matrix bestMSoFar; GetAllIndependent (bestMSoFar); hyFloat prec = Minimum (diffs[0], diffs[1]), - grad_precision = diffs[0] + diffs[1]; + grad_precision = Maximum(diffs[0],diffs[1]); prec = Minimum (Maximum (prec, precision), 1.); @@ -5182,7 +5328,7 @@ _Matrix* _LikelihoodFunction::Optimize (_AssociativeList const * options) for (unsigned long i=0UL; ivoid { for (long i=0L; i 100) { @@ -6314,8 +6462,9 @@ void _LikelihoodFunction::ComputeGradient (_Matrix& gradient, hyFloat& gradi hyFloat currentValue = GetIthIndependent(index), ub = GetIthIndependentBound(index,false)-currentValue, lb = currentValue-GetIthIndependentBound(index,true), - testStep = currentValue * gradientStep;//MAX(currentValue * gradientStep,gradientStep); - + testStep = MAX(currentValue * gradientStep,gradientStep); + + //printf ("%ld %s %20.18g\n", index, GetIthIndependentName (index)->get_str(), currentValue); hyFloat check_vv = currentValue; @@ -6332,11 +6481,11 @@ void _LikelihoodFunction::ComputeGradient (_Matrix& gradient, hyFloat& gradi } } - + if (!CheckEqual (testStep,0.0)) { - /*if (verbosity_level > 100) { - printf ("Gradient step for %s is %.16g @%.16g %\n", GetIthIndependentVar(index)->GetName()->get_str(), testStep, currentValue); - }*/ + //if (verbosity_level > 100) { + //printf ("Gradient step for %s is %.16g @%.16g %\n", GetIthIndependentVar(index)->GetName()->get_str(), testStep, currentValue); + //} SetIthIndependent(index,currentValue+testStep); hyFloat dF = Compute(); gradient[index]=(dF-funcValue)/testStep;// * DerivativeCorrection (index, currentValue); @@ -6349,9 +6498,9 @@ void _LikelihoodFunction::ComputeGradient (_Matrix& gradient, hyFloat& gradi printf ("Negative value stashed %15.12lg\n", currentValue); } hyFloat check_vv = GetIthIndependent(index);*/ - if (verbosity_level > 150) { - printf ("_LikelihoodFunction::ComputeGradient %ld\t%s\t @%20.18g\t f(x) = %20.18g, f(x+h) = %20.18g, h = %20.18g, correction = %20.18g, grad = %20.18g\n", index, GetIthIndependentName(index)->get_str(), currentValue, funcValue, dF, testStep, DerivativeCorrection (index, currentValue), gradient.theData[index]); - } + //if (verbosity_level > 150) { + //printf ("_LikelihoodFunction::ComputeGradient %ld\t%s\t @%20.18g\t f(x) = %20.18g, f(x+h) = %20.18g, delta(F) = %20.18g, h = %20.18g, correction = %20.18g, grad = %20.18g\n", index, GetIthIndependentName(index)->get_str(), currentValue, funcValue, dF, dF-funcValue, testStep, DerivativeCorrection (index, currentValue), gradient.theData[index]); + //} SetIthIndependent(index,currentValue); } else { gradient[index]= 0.; @@ -6542,8 +6691,8 @@ hyFloat _LikelihoodFunction::ConjugateGradientDescent (hyFloat step_precision if (gradL > 0.0) { - current_direction = gradient * (1./gradL); - gradient = current_direction; + current_direction = gradient;// * (1./gradL); + //gradient = current_direction; for (long index = 0; index< MAX (dim, 10) && index < iterationLimit; index++, currentPrecision*=0.25) { @@ -6570,33 +6719,30 @@ hyFloat _LikelihoodFunction::ConjugateGradientDescent (hyFloat step_precision ComputeGradient (gradient, gradientStep, bestVal, freeze, 1, false); gradL = gradient.AbsValue (); - if (CheckEqual(gradL,0.0)) { + /*if (CheckEqual(gradL,0.0)) { break; } else { gradient *= 1./gradL; - } + }*/ previous_direction = current_direction; hyFloat beta = 0., scalar_product = 0.; // use Polak–Ribière direction - for (unsigned long i = 0UL; i < dim; i++) { + /*for (unsigned long i = 0UL; i < dim; i++) { scalar_product += previous_gradient.theData[i] * previous_gradient.theData[i]; beta += gradient.theData[i] * ( previous_gradient.theData[i] - gradient.theData[i]); - } + }*/ - LoggerAddGradientPhase (line_search_precision, beta, scalar_product); - LoggerAllVariables (); - LoggerLogL (maxSoFar); // use Dai–Yuan - /* + for (unsigned long i = 0UL; i < dim; i++) { beta += gradient.theData[i] * gradient.theData[i]; scalar_product += previous_direction.theData[i] * ( gradient.theData[i] - previous_gradient.theData[i]); } beta = -beta; - */ + // Hestenes-Stiefel /* @@ -6607,6 +6753,10 @@ hyFloat _LikelihoodFunction::ConjugateGradientDescent (hyFloat step_precision beta = -beta; */ + LoggerAddGradientPhase (line_search_precision, beta, scalar_product); + LoggerAllVariables (); + LoggerLogL (maxSoFar); + beta /= scalar_product; beta = MAX (beta, 0.0); previous_direction = current_direction; diff --git a/src/core/matrix.cpp b/src/core/matrix.cpp index 5fa1a3730..3c4f620b9 100644 --- a/src/core/matrix.cpp +++ b/src/core/matrix.cpp @@ -2501,25 +2501,6 @@ bool _Matrix::ProcessFormulas (long& stackLength, _AVLList& varList, _S _Formula ** theFormulas = (_Formula**)theData; bool isGood = true; - - auto process_formula = [&] (_Formula * this_formula) -> bool { - if (runAll || this_formula->AmISimple(stackLength,varList)) { - _String * flaString = (_String*)this_formula->toStr(kFormulaStringConversionNormal, nil,true); - long fref = flaStrings.Insert(flaString,newFormulas.lLength); - if (fref < 0) { - references << flaStrings.GetXtra (-fref-1); - DeleteObject (flaString); - } else { - newFormulas << (long)this_formula; - references << fref; - } - - } else { - isGood = false; - return false; - } - return true; - }; if (theIndex) { for (long i = 0L; iProcessFormulaForConverstionToSimple (stackLength, varList, newFormulas, references, flaStrings, runAll)) { + isGood = false; break; } } else { @@ -2546,7 +2528,8 @@ bool _Matrix::ProcessFormulas (long& stackLength, _AVLList& varList, _S references << -1; continue; } - if (! process_formula (thisFormula)) { + if (! thisFormula->ProcessFormulaForConverstionToSimple (stackLength, varList, newFormulas, references, flaStrings, runAll)) { + isGood = false; break; } @@ -2718,8 +2701,16 @@ void _Matrix::MakeMeSimple (void) { memcpy (cmd->formulaRefs, references.list_data, allocation_size); cmd->formulaValues = new hyFloat [newFormulas.lLength]; cmd->formulasToEval.Duplicate (&newFormulas); +/* + printf ( +"\nConverted a matrix to simple expressions. \ +\n\tStack depth: %ld\ +\n\tIndependent variables: %ld\ +\n\tUnique expressions: %ld\n", stackLength, varList.countitems(), cmd->formulasToEval.countitems()); +*/ } + } } //_____________________________________________________________________________________________ diff --git a/src/core/topology.cpp b/src/core/topology.cpp index ec6bca60b..f871d92d5 100644 --- a/src/core/topology.cpp +++ b/src/core/topology.cpp @@ -259,12 +259,14 @@ struct _any_char_in_set { const _TreeTopologyParseSettings _TreeTopology::CollectParseSettings (void) { _String const kInternalNodePrefix("INTERNAL_NODE_PREFIX"), - kIgnoreUserINames ("IGNORE_INTERNAL_NODE_LABELS"); + kIgnoreUserINames ("IGNORE_INTERNAL_NODE_LABELS"), + kDoNotAutoRename ("TREE_DO_NOT_AUTO_RENAME"); _TreeTopologyParseSettings parse_settings; parse_settings.auto_convert_lengths = EnvVariableTrue(automatically_convert_branch_lengths); parse_settings.accept_user_lengths = EnvVariableTrue(accept_branch_lengths); + parse_settings.auto_rename_nodes = !EnvVariableTrue(kDoNotAutoRename); parse_settings.ingore_user_inode_names = EnvVariableTrue(kIgnoreUserINames); HBLObjectRef user_node_name = EnvVariableGet(kInternalNodePrefix, STRING); @@ -400,7 +402,10 @@ _AssociativeList* _TreeTopology::MainTreeConstructor (_String const& parms, if (mapped_name) { //printf ("%s => %s\n", nodeName.get_str(), mapped_name->get_str().get_str()); nodeName = _String (mapped_name->get_str()); - } + } /*else { + //ObjectToConsole(mapping); + printf ("%s => NOTHING\n", nodeName.get_str()); + }*/ } // handle bootstrap support for this node diff --git a/src/core/tree.cpp b/src/core/tree.cpp index 56504acaa..f407fc50a 100644 --- a/src/core/tree.cpp +++ b/src/core/tree.cpp @@ -505,9 +505,13 @@ _String _TheTree::FinalizeNode (node* nodie, long number , _String node nodeName = settings.inode_prefix & number; } else { if (!nodeName.IsValidIdentifier(fIDAllowFirstNumeric)) { - _String new_name = nodeName.ConvertToAnIdent(fIDAllowFirstNumeric); - ReportWarning (_String ("Automatically renamed ") & nodeName.Enquote() & " to " & new_name.Enquote() & " in order to create a valid HyPhy identifier"); - nodeName = new_name; + if (settings.auto_rename_nodes) { + _String new_name = nodeName.ConvertToAnIdent(fIDAllowFirstNumeric); + ReportWarning (_String ("Automatically renamed ") & nodeName.Enquote() & " to " & new_name.Enquote() & " in order to create a valid HyPhy identifier"); + nodeName = new_name; + } else { + HandleApplicationError(nodeName.Enquote() & " is not a valid HyPhy identifier (auto-remaing not allowed)."); + } } } diff --git a/src/core/tree_evaluator.cpp b/src/core/tree_evaluator.cpp index 4f78f10ed..358e4bd3c 100644 --- a/src/core/tree_evaluator.cpp +++ b/src/core/tree_evaluator.cpp @@ -157,9 +157,9 @@ template inline bool __ll_handle_conditional_array_initialization ( long long siteState; if (setBranch != nodeCode + iNodes) { siteState = lNodeFlags[nodeCode*siteCount + siteOrdering.list_data[siteID]] ; - /*if (likeFuncEvalCallCount == 328 && siteID == 30) { - fprintf (stderr, "Site-state @ %ld: %d\n", nodeCode, siteState); - }*/ + //if (siteID == 524) { + // fprintf (stderr, "Site-state @ %ld: %d\n", nodeCode, siteState); + //} } else { siteState = setBranchTo[siteOrdering.list_data[siteID]] ; } @@ -190,6 +190,14 @@ template inline bool __ll_handle_conditional_array_initialization ( long for (long k = ub; k < D; k++) { parentConditionals[k] *= tMatrix[siteState+D*k]; } + + /* + if (parentConditionals[siteState] == 0.) { + if (siteID == 524) { + fprintf (stderr, "%ld/%ld / %g\n", siteState, siteID, parentConditionals[siteState]); + } + } + */ #else #pragma GCC unroll 4 for (long k = 0L; k < D; k++) { @@ -3052,8 +3060,8 @@ hyFloat _TheTree::ComputeTreeBlockByBranch ( _SimpleList //if (setBranch >=0 ) // printf ("\nSet to %d (%s)\n", setBranch, ((_CalcNode*) flatTree (setBranch))->GetName()->get_str()); - /*if (likeFuncEvalCallCount == 328) { - fprintf (stderr, "\nSite ID: %ld\n%s\n", siteOrdering.get(30), theFilter->GetColumn(siteOrdering.get(30)));; + /*if (likeFuncEvalCallCount == 0) { + fprintf (stderr, "\nSite ID: %ld\n%s\n%s\n%s\n", siteOrdering.get(0), theFilter->GetColumn(siteOrdering.get(0)), theFilter->GetColumn(siteOrdering.get(0)+1), theFilter->GetColumn(siteOrdering.get(0)+2)); }*/ for (unsigned long nodeID = 0; nodeID < updateNodes.lLength; nodeID++) { @@ -3295,7 +3303,7 @@ hyFloat _TheTree::ComputeTreeBlockByBranch ( _SimpleList /* #pragma omp critical - if (siteID == siteFrom && nodeID == 25) { + if (siteID == 524) { printf ("(%ld/%ld)%s\n", likeFuncEvalCallCount, siteFrom, currentTreeNode->GetName()->get_str()); } */ @@ -3310,15 +3318,18 @@ hyFloat _TheTree::ComputeTreeBlockByBranch ( _SimpleList } -//#pragma omp critical - /*if (1 || hy_env::EnvVariableTrue("UBER_VERBOSE_DEBUG")) { - if (siteOrdering.list_data[siteID] == 2 && nodeID == 75) { +/* +#pragma omp critical + if (1 || hy_env::EnvVariableTrue("UBER_VERBOSE_DEBUG")) { + if (siteID == 524) { + //if (siteOrdering.list_data[siteID] == 2 && nodeID == 75) { //printf ("%ld\t%ld\t%ld\t%ld\t%16.12g\n",nodeID, siteFrom, siteTo, siteID, sum); for (int kk = 0; kk < 61; kk++) { printf ("%s (%ld)\t%ld\t%ld\t%d\t%16.12g\t%16.12g\t%16.12g\n", currentTreeNode->GetName()->get_str(),nodeID, siteFrom, siteID, kk, parentConditionals[kk], childVector[kk], mvs[kk]); } } - }*/ + } +*/ #ifdef _SLKP_USE_APPLE_BLAS_2 cblas_dgemv(CblasRowMajor, @@ -3338,16 +3349,18 @@ hyFloat _TheTree::ComputeTreeBlockByBranch ( _SimpleList #endif sum += _hy_vvmult_sum<61> (parentConditionals, mvs); - -//#pragma omp critical - /*if (likeFuncEvalCallCount == 1 || hy_env::EnvVariableTrue("UBER_VERBOSE_DEBUG")) { - if (siteOrdering.list_data[siteID] == 2 && nodeID == 75) { +/* +#pragma omp critical + if (1 || hy_env::EnvVariableTrue("UBER_VERBOSE_DEBUG")) { + if (siteID == 524) { + //if (siteOrdering.list_data[siteID] == 2 && nodeID == 75) { //printf ("%ld\t%ld\t%ld\t%ld\t%16.12g\n",nodeID, siteFrom, siteTo, siteID, sum); for (int kk = 0; kk < 61; kk++) { - printf ("%s\t%ld\t%ld\t%d\t%16.12g\t%16.12g\t%16.12g\n", currentTreeNode->GetName()->get_str(), siteFrom, siteID, kk, parentConditionals[kk], childVector[kk], mvs[kk]); + printf ("%s (%d)\t%ld\t%ld\t%d\t%16.12g\t%16.12g\t%16.12g\n", currentTreeNode->GetName()->get_str(), nodeID, siteFrom, siteID, kk, parentConditionals[kk], childVector[kk], mvs[kk]); } } - }*/ + } +*/ @@ -3498,7 +3511,7 @@ hyFloat _TheTree::ComputeTreeBlockByBranch ( _SimpleList if (accumulator <= 0.0) #pragma omp critical { - //printf ("BAILING WITH INFINITY %ld\n", siteID); + //printf ("BAILING WITH INFINITY %ld / %ld\n", siteID, siteOrdering.list_data[siteID]); hy_global::ReportWarning (_String("Site ") & (1L+siteOrdering.list_data[siteID]) & " evaluated to a 0 probability in ComputeTreeBlockByBranch with category " & catID); } } else { @@ -3506,7 +3519,7 @@ hyFloat _TheTree::ComputeTreeBlockByBranch ( _SimpleList result = -INFINITY; #pragma omp critical { - //printf ("BAILING WITH INFINITY %ld\n", siteID); + //printf ("BAILING WITH INFINITY %ld / %ld (%ld)\n", siteID, siteOrdering.list_data[siteID], siteOrdering.lLength); hy_global::ReportWarning (_String("Site ") & (1L+siteOrdering.list_data[siteID]) & " evaluated to a 0 probability in ComputeTreeBlockByBranch"); } break; diff --git a/src/core/variable.cpp b/src/core/variable.cpp index e1e2285a0..f67dbaa0c 100644 --- a/src/core/variable.cpp +++ b/src/core/variable.cpp @@ -264,9 +264,14 @@ HBLObjectRef _Variable::Compute (void) { auto update_var_value = [this] () -> void { if (!varValue || varFormula->HasChanged()) { - HBLObjectRef new_value = (HBLObjectRef)varFormula->Compute()->makeDynamic(); + HBLObjectRef new_value = (HBLObjectRef)varFormula->Compute(); + if (varValue && new_value->ObjectClass () == NUMBER) { + if (varValue->ObjectClass () == NUMBER && varValue->CanFreeMe()) { + ((_Constant*)varValue)->SetValue(((_Constant*)new_value)->Value()); + } + } DeleteObject (varValue); - varValue = new_value; + varValue = (HBLObjectRef)new_value->makeDynamic(); } }; @@ -509,6 +514,7 @@ void _Variable::SetNumericValue (hyFloat v) // set the value of the var to a nu varFlags &= HY_VARIABLE_SET; varFlags |= HY_VARIABLE_CHANGED; theValue = v; + if (theValueupperBound) { @@ -706,7 +712,12 @@ void _Variable::SetFormula (_Formula& theF) { varFormula->Duplicate (right_hand_side); // mod 20060125 added a call to simplify constants + + //_String * fs = (_String*)varFormula->toStr(kFormulaStringConversionNormal); + //NLToConsole(); StringToConsole(*fs); DeleteObject (fs); varFormula->SimplifyConstants (); + //fs = (_String*)varFormula->toStr(kFormulaStringConversionNormal); + //BufferToConsole (" : "); StringToConsole(*fs); DeleteObject (fs); // also update the fact that this variable is no longer independent in all declared // variable containers which hold references to this variable diff --git a/tests/hbltests/libv3/BUSTED.wbf b/tests/hbltests/libv3/BUSTED.wbf index 1a44ff4e1..d42f6e2a6 100644 --- a/tests/hbltests/libv3/BUSTED.wbf +++ b/tests/hbltests/libv3/BUSTED.wbf @@ -1,7 +1,7 @@ GetString (version, HYPHY_VERSION, 0); if (+version >= 2.4) { - LoadFunctionLibrary ("SelectionAnalyses/BUSTED.bf", {"--code" : "Universal", "--alignment" : PATH_TO_CURRENT_BF + "data/CD2.nex", "--branches" : "GROUP1", "--srv" : "No", "--starting-points" : "5"}); + LoadFunctionLibrary ("SelectionAnalyses/BUSTED.bf", {"--code" : "Universal", "--alignment" : PATH_TO_CURRENT_BF + "data/CD2.nex", "--branches" : "GROUP1", "--srv" : "Yes", "--starting-points" : "5"}); } else { LoadFunctionLibrary ("SelectionAnalyses/BUSTED.bf", {"0" : "Universal", "1" : PATH_TO_CURRENT_BF + "data/CD2.nex", "2" : "GROUP1", "3" : "Yes", "4" : "0.1"}); @@ -15,14 +15,14 @@ LoadFunctionLibrary ("libv3/IOFunctions.bf"); assert (check_value ( - ((busted.json["fits"])["Unconstrained model"])["Log Likelihood"], -3435.55, 0.001), "Incorrect log-likelihood for the full model"); + ((busted.json["fits"])["Unconstrained model"])["Log Likelihood"], -3413.83, 0.001), "Incorrect log-likelihood for the full model"); assert (check_value ( - ((busted.json["test results"])["p-value"]),0.18531, 0.01), "p-value for the test is incorrect"); + ((busted.json["test results"])["p-value"]),0.2264, 0.01), "p-value for the test is incorrect"); assert (check_value ( - +(busted.json["Evidence Ratios"])["optimized null"],189.46, 0.01), "Incorrect sum of evidence ratios"); + +(busted.json["Evidence Ratios"])["optimized null"],188.49, 0.01), "Incorrect sum of evidence ratios");