diff --git a/res/TemplateBatchFiles/SelectionAnalyses/contrast-fel.bf b/res/TemplateBatchFiles/SelectionAnalyses/contrast-fel.bf index d56507664..152f7f866 100644 --- a/res/TemplateBatchFiles/SelectionAnalyses/contrast-fel.bf +++ b/res/TemplateBatchFiles/SelectionAnalyses/contrast-fel.bf @@ -396,12 +396,26 @@ fel.site_patterns = alignments.Extract_site_patterns ((fel.filter_specification[ // alpha = alpha_scaler * branch_length // beta = beta_scaler_test * branch_length or beta_nuisance_test * branch_length -utility.ForEach (fel.case_respecting_node_names, "_node_", - '_node_class_ = (fel.selected_branches[fel.partition_index])[_node_]; - _beta_scaler = fel.scaler_parameter_names[_node_class_]; - fel.apply_proportional_site_constraint ("fel.site_tree", _node_, fel.alpha, fel.beta, fel.alpha.scaler, _beta_scaler, (( fel.final_partitioned_mg_results[terms.branch_length])[fel.partition_index])[_node_]); - '); +fel.lengths_by_class = {}; + +for (_node_; in; fel.case_respecting_node_names) { + _node_class_ = (fel.selected_branches[fel.partition_index])[_node_]; + _beta_scaler = fel.scaler_parameter_names[_node_class_]; + fel.lengths_by_class[_node_class_] += fel.apply_proportional_site_constraint ("fel.site_tree", _node_, fel.alpha, fel.beta, fel.alpha.scaler, _beta_scaler, (( fel.final_partitioned_mg_results[terms.branch_length])[fel.partition_index])[_node_]); +} +fel.ignorable = {}; +for (fel.t; in; fel.branches.testable) { + if ( fel.lengths_by_class[fel.t] == 0) { + fprintf(stdout, "\n-------\n", io.FormatLongStringToWidth( + ">[WARNING] The cumulative branch length in the _" + fel.t + "_ class is 0. + Rates along these branches are not identifiable; testing will not be formally conducted (all p-values set to 1).", 72), + "\n-------\n"); + fel.ignorable [fel.scaler_parameter_names[fel.t]] = 1; + } else { + fel.ignorable [fel.scaler_parameter_names[fel.t]] = 0; + } +} // create the likelihood function for this site @@ -424,7 +438,7 @@ fel.queue = mpi.CreateQueue ({"LikelihoodFunctions": {{"fel.site_likelihood"}}, "Models" : {{"fel.site.mg_rev"}}, "Headers" : {{"libv3/all-terms.bf","libv3/tasks/ancestral.bf", "libv3/convenience/math.bf"}}, "Functions" : {{"fel.apply_proportional_site_constraint"}}, - "Variables" : {{"terms.fel.test_keys","fel.permutations", "fel.alpha","fel.beta","fel.alpha.scaler","terms.fel.permutation","fel.final_partitioned_mg_results","fel.srv","fel.site_tested_classes","fel.scaler_parameter_names","fel.branches.testable","fel.branches.has_background","fel.alpha.scaler","terms.fel.pairwise","fel.branch_class_counter","fel.report.test_count", "fel.p_value","fel.site.permutations"}} + "Variables" : {{"terms.fel.test_keys","fel.permutations","fel.ignorable","fel.alpha","fel.beta","fel.alpha.scaler","terms.fel.permutation","fel.final_partitioned_mg_results","fel.srv","fel.site_tested_classes","fel.scaler_parameter_names","fel.branches.testable","fel.branches.has_background","fel.alpha.scaler","terms.fel.pairwise","fel.branch_class_counter","fel.report.test_count", "fel.p_value","fel.site.permutations"}} }); @@ -644,6 +658,8 @@ function fel.apply_proportional_site_constraint (tree_name, node_name, alpha_par `node_name`.`alpha_parameter` := (`alpha_factor`) * fel.branch_length__; `node_name`.`beta_parameter` := (`beta_factor`) * fel.branch_length__; "); + + return fel.branch_length; } //---------------------------------------------------------------------------------------- @@ -724,8 +740,6 @@ lfunction fel.handle_a_site (lf, filter_data, partition_index, pattern_info, mod '); } - - snapshot = estimators.TakeLFStateSnapshot (lf); alternative = estimators.ExtractMLEs (lf, model_mapping); alternative [utility.getGlobalValue("terms.fit.log_likelihood")] = results[1][0]; @@ -757,8 +771,8 @@ lfunction fel.handle_a_site (lf, filter_data, partition_index, pattern_info, mod } else { for (_gname_; in; ^"fel.branches.testable") { _pname_ = (^"fel.scaler_parameter_names")[_gname_]; - if (_pname_ != ref_parameter) { - ^_pname_ := ^ref_parameter; + if (_pname_ != ref_parameter && (^"fel.ignorable")[_pname_] == 0) { + parameters.SetConstraint (_pname_,ref_parameter, ""); } } } @@ -773,12 +787,17 @@ lfunction fel.handle_a_site (lf, filter_data, partition_index, pattern_info, mod for (v2 = v + 1; v2 < testable; v2+=1) { v1n = (^"fel.branches.testable")[v]; v2n = (^"fel.branches.testable")[v2]; - estimators.RestoreLFStateFromSnapshot (lf_id, snapshot); - parameters.SetConstraint ((^"fel.scaler_parameter_names")[v1n],(^"fel.scaler_parameter_names")[v2n], ""); - Optimize (results, ^lf, {"OPTIMIZATION_METHOD" : "nedler-mead"}); - pairwise[v1n + "|" + v2n] = estimators.ExtractMLEs (lf, model_mapping); - (pairwise[v1n + "|" + v2n])[utility.getGlobalValue("terms.fit.log_likelihood")] = results[1][0]; + + if ((^"fel.ignorable")[(^"fel.scaler_parameter_names")[v1n]] || (^"fel.ignorable")[(^"fel.scaler_parameter_names")[v2n]]) { + //console.log (v1n + "|" + v2n + " is ignorable"); + (pairwise[v1n + "|" + v2n]) = alternative; + } else { + parameters.SetConstraint ((^"fel.scaler_parameter_names")[v1n],(^"fel.scaler_parameter_names")[v2n], ""); + Optimize (results, ^lf, {"OPTIMIZATION_METHOD" : "nedler-mead"}); + pairwise[v1n + "|" + v2n] = estimators.ExtractMLEs (lf, model_mapping); + (pairwise[v1n + "|" + v2n])[utility.getGlobalValue("terms.fit.log_likelihood")] = results[1][0]; + } } } } else { diff --git a/src/core/batchlanhelpers.cpp b/src/core/batchlanhelpers.cpp index d2a886573..b8bfa14d1 100644 --- a/src/core/batchlanhelpers.cpp +++ b/src/core/batchlanhelpers.cpp @@ -68,14 +68,18 @@ _String const _HYGenerateANameSpace () { //____________________________________________________________________________________ void _HYClearANameSpace (const _String& nm) { + /*BufferToConsole("Deleting key: "); + BufferToConsole(nm.get_str()); + NLToConsole(); + ObjectToConsole(&_HY_HBL_Namespaces); + NLToConsole();*/ _HY_HBL_Namespaces.Delete(nm); } //____________________________________________________________________________________ -void ReadModelList(void) -{ +void ReadModelList(void) { if (templateModelList.empty() == false) return; _String modelListFile (GetStandardDirectory (HY_HBL_DIRECTORY_TEMPLATE_MODELS) & "models.lst"), diff --git a/src/core/global_things.cpp b/src/core/global_things.cpp index 869ed8556..7749410de 100644 --- a/src/core/global_things.cpp +++ b/src/core/global_things.cpp @@ -121,7 +121,7 @@ namespace hy_global { kErrorStringDatasetRefIndexError ("Dataset index reference out of range"), 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"), - kHyPhyVersion = _String ("2.5.22"), + kHyPhyVersion = _String ("2.5.23"), kNoneToken = "None", kNullToken = "null", diff --git a/src/core/likefunc.cpp b/src/core/likefunc.cpp index 627b945b1..e307b6174 100644 --- a/src/core/likefunc.cpp +++ b/src/core/likefunc.cpp @@ -4101,6 +4101,9 @@ _Matrix* _LikelihoodFunction::Optimize (_AssociativeList const * options) #endif + _variables_changed_during_last_compute = new _SimpleList (); + variables_changed_during_last_compute = new _AVLList (_variables_changed_during_last_compute); + #ifdef __HYPHYMPI__ if (hy_mpi_node_rank == 0) { #endif @@ -4109,9 +4112,6 @@ _Matrix* _LikelihoodFunction::Optimize (_AssociativeList const * options) } #endif - _variables_changed_during_last_compute = new _SimpleList (); - variables_changed_during_last_compute = new _AVLList (_variables_changed_during_last_compute); - bool skipCG = ! CheckEqual (get_optimization_setting (kSkipConjugateGradient, 0.0), 0.0), keepStartingPoint = ! CheckEqual (get_optimization_setting (kUseLastResults, 0.0), 0.0), diff --git a/src/core/trie.cpp b/src/core/trie.cpp index 8dba0da31..b64b52a23 100644 --- a/src/core/trie.cpp +++ b/src/core/trie.cpp @@ -344,10 +344,11 @@ bool _Trie::Delete (const _String& key){ emptySlots << history.list_data[k]; payload.list_data[history.list_data[k]] = 0L; parents.list_data[history.list_data[k]] = -1L; - _SimpleList * parentList = ((_SimpleList**)list_data)[history.list_data[k-1]]; + _SimpleList * parentList = ((_SimpleList**)list_data)[k ? history.list_data[k-1] : 0]; unsigned long parentNode = parentList->FindStepping (history.list_data[k],2, 1) - 1; parentList->Delete (parentNode); parentList->Delete (parentNode); + DeleteObject (current_list); ((_SimpleList**)list_data)[history.list_data[k]] = nil; inserted_values--; @@ -450,7 +451,9 @@ BaseRef _Trie::toStr(unsigned long) { } } else { _String * this_string = RetrieveStringFromPath(traversal_history, &alph); - (*serialized) << '"' << this_string << "\":" << _String (GetValue (traversal_history.list_data[traversal_history.lLength-2])); + (*serialized) << '"'; + (*serialized).AppendNewInstance(this_string); + (*serialized) << "\":" << _String (GetValue (traversal_history.list_data[traversal_history.lLength-2])); if (doComma) { (*serialized) << ','; } else { diff --git a/src/mains/unix.cpp b/src/mains/unix.cpp index 13cbb4d7a..020c7d6fe 100644 --- a/src/mains/unix.cpp +++ b/src/mains/unix.cpp @@ -694,51 +694,6 @@ void SetStatusLineUser (_String const s) { #ifndef __UNITTEST__ int main (int argc, char* argv[]) { - /*long N = 1000L; - long repeats = 20L; - long density = 0.25 * RAND_MAX; - - auto start = std::chrono::high_resolution_clock::now(); - - srand(time(0)); - - _Matrix::switchThreshold = atoi (argv[1]); - - printf ("%ld\n", _Matrix::switchThreshold); - - for (long i = 0; i < N; i++) { - _Matrix test (61,61,false,true); - test.ForEachCellNumeric([density] (hyFloat& e, long i, long r, long c) -> void { - if (r != c && rand() < density) { - e = genrand_real1 (); - } else { - e = 0.0; - } - }); - for (long r = 0; r < 61; r++) { - double s = 0.; - for (long c = 0; c < 61; c++) { - s += test (r, c); - } - test.theData[r*61+r] = -s; - } - - test.AmISparse(); - //ObjectToConsole(&test); - - //printf ("%d\n", test.is_dense()); - for (long j = 0; j < repeats; j++) { - DeleteObject (test.Exponentiate()); - //return 0; - } - } - - auto end = std::chrono::high_resolution_clock::now(); - std::chrono::duration elapsed = end - start; - printf ("%15.12g s\n", elapsed.count()); - - return 0;*/ - #ifdef _COMPARATIVE_LF_DEBUG_DUMP FILE * comparative_lf_debug_matrix_content_file = doFileOpen (_COMPARATIVE_LF_DEBUG_DUMP, "w");