Merge pull request #1926 from veg/develop

2.5.86

Author: spond

res/TemplateBatchFiles/SelectionAnalyses/RELAX.bf

@@ -25,6 +25,8 @@ LoadFunctionLibrary("libv3/models/codon/BS_REL.bf");
 LoadFunctionLibrary("libv3/convenience/math.bf");
 LoadFunctionLibrary("libv3/tasks/mpi.bf");
 LoadFunctionLibrary("libv3/models/rate_variation.bf");
+LoadFunctionLibrary("libv3/tasks/ancestral.bf");
+
 
 
 
@@ -179,7 +181,14 @@ namespace relax {
     LoadFunctionLibrary ("modules/grid_compute.ibf");
 }
 
-KeywordArgument ("output", "Write the resulting JSON to this file (default is to save to the same path as the alignment file + 'RELAX.json')", relax.codon_data_info [terms.json.json]);
+
+relax.default_output = relax.codon_data_info [terms.json.json];
+
+if (relax.multiple_files) {
+    relax.default_output = (relax.codon_data_info [relax.file_list[0]])[terms.json.json];
+}
+
+KeywordArgument ("output", "Write the resulting JSON to this file (default is to save to the same path as the alignment file + 'RELAX.json')", relax.default_output);
 
 relax.codon_data_info [terms.json.json] = io.PromptUserForFilePath ("Save the resulting JSON file to");
 
@@ -1180,6 +1189,7 @@ function relax.FitMainTestPair (prompt) {
                 relax.stashLF = estimators.TakeLFStateSnapshot (relax.alternative_model.fit[terms.likelihood_function]);
 
 			} else {
+			    io.ReportProgressMessageMD("RELAX", "alt-2", "> Restoring the original set of parameters (alternative optimization did not yield an improvement).");
 			    estimators.RestoreLFStateFromSnapshot(relax.alternative_model.fit[terms.likelihood_function], relax.take1_snapshot);
 			}
 
@@ -1218,6 +1228,16 @@ function relax.FitMainTestPair (prompt) {
 								relax.distribution_for_json,
 								relax.display_orders[relax.alternative_name]
 							);
+							
+    relax.json [^"terms.substitutions"] = {};
+    
+    for (_partition_, _selection_; in; relax.selected_branches) {
+        relax.ancestral_info = ancestral.build (relax.alternative_model.fit[terms.likelihood_function],0+_partition_,FALSE);
+        (relax.json [^"terms.substitutions"] )[_partition_] = ancestral.ComputeCompressedSubstitutions (relax.ancestral_info);
+        DeleteObject (relax.ancestral_info);
+    }
+
+
 
 	selection.io.json_store_branch_attribute(relax.json, relax.alternative_name, terms.branch_length, relax.display_orders[relax.alternative_name],
 												 0,
@@ -1231,9 +1251,15 @@ function relax.FitMainTestPair (prompt) {
     relax.tree_ids = estimators.LFObjectGetTrees (relax.alternative_model.fit[terms.likelihood_function]);
 
     estimators.RestoreLFStateFromSnapshot (relax.alternative_model.fit[terms.likelihood_function],relax.stashLF);
+    
 
     utility.ToggleEnvVariable ("KEEP_OPTIMAL_ORDER", TRUE);
         relax.er_report.tagged_sites    = {};
+        relax.json [terms.json.site_logl] = {};
+        relax.json [terms.json.evidence_ratios]  = {};
+        (relax.json [terms.json.site_logl])[terms.json.unconstrained] = selection.ComputeSiteLikelihoods (relax.alternative_model.fit[terms.likelihood_function]);
+
+
         for (_partition_, _selection_; in; relax.selected_branches) {
 
 
@@ -1254,10 +1280,8 @@ function relax.FitMainTestPair (prompt) {
             relax.ll_by_branch = {};
             relax.model_by_branch = {};
 
-            
+            (relax.json [terms.json.site_logl])[_partition_] = {};
             relax.altLL = relax.alternative_model.fit [terms.fit.log_likelihood];
-            relax.json [terms.json.site_logl] = {};
-            (relax.json [terms.json.site_logl])[terms.json.unconstrained] = selection.ComputeSiteLikelihoods (relax.alternative_model.fit[terms.likelihood_function]);
 
             for (_b_,_ignore_; in;  relax.tested_branches) {
                  GetString    (relax._current_model, ^_b_, -5);                 
@@ -1272,22 +1296,21 @@ function relax.FitMainTestPair (prompt) {
                  SetParameter ( ^_b_ , MODEL, ^relax._current_model);
             }
 
-
-            relax.json [terms.json.evidence_ratios] = {terms.json.by_branch : relax.ll_by_branch};
+            (relax.json [terms.json.evidence_ratios])[_partition_] = {terms.json.by_branch : relax.ll_by_branch};
 
             for (_b_,_ignore_; in;  relax.tested_branches) {
                  SetParameter ( ^_b_ , MODEL, ^relax.reference_model_namespace);
             }
             LFCompute (^(relax.alternative_model.fit[terms.likelihood_function]),LF_START_COMPUTE);
             LFCompute (^(relax.alternative_model.fit[terms.likelihood_function]),relax.ll);
             LFCompute (^(relax.alternative_model.fit[terms.likelihood_function]),LF_DONE_COMPUTE);
-            (relax.json [terms.json.site_logl])[terms.json.constrained] = selection.ComputeSiteLikelihoods (relax.alternative_model.fit[terms.likelihood_function]);
+            ((relax.json [terms.json.site_logl])[_partition_])[terms.json.constrained] = selection.ComputeSiteLikelihoods (relax.alternative_model.fit[terms.likelihood_function]);
 
             for (_b_,_ignore_; in;  relax.tested_branches) {
                  SetParameter ( ^_b_ , MODEL, ^(relax.model_by_branch[_b_]));
             }
 
-            (relax.json [terms.json.evidence_ratios])[terms.json.constrained] = selection.EvidenceRatios ( (relax.json [terms.json.site_logl])[terms.json.unconstrained],  (relax.json [terms.json.site_logl])[terms.json.constrained]);
+            ((relax.json [terms.json.site_logl])[_partition_])[terms.json.constrained] = selection.EvidenceRatios ((relax.json [terms.json.site_logl])[terms.json.unconstrained],  ((relax.json [terms.json.site_logl])[_partition_])[terms.json.constrained]);
 
 
 

res/TemplateBatchFiles/libv3/tasks/genetic_code.bf

@@ -91,6 +91,21 @@ lfunction genetic_code.ComputeCodonCodeToStringMap(genCode) {
 
 /*----------------------------------------------------------------------------------------------------------*/
 
+lfunction genetic_code.ComputeStringToInteger (genCode) {
+    _codonMap = {};
+    _nucLetters = ^"genetic_code.nucleotides";
+    _cidx = 0;
+    for (_idx = 0; _idx < Columns(genCode); _idx += 1) {
+        if (genCode[_idx] != ^ "genetic_code.stop_code") {
+            _codonMap[(_nucLetters[_idx$16] + _nucLetters[(_idx % 16) $4] + _nucLetters[_idx % 4])] = _cidx;
+            _cidx += 1;
+        }
+    }
+    return _codonMap;
+}
+
+/*----------------------------------------------------------------------------------------------------------*/
+
 lfunction genetic_code.ComputeCodonCodeToStringMapStop (genCode) {
 	_codonMap = {};
 	_nucLetters =  ^"genetic_code.nucleotides";

res/TemplateBatchFiles/libv3/tasks/trees.bf

@@ -660,7 +660,11 @@ lfunction trees.ParentMap (topology_object) {
     for (i,k; in; parent_nodes) {
         if (node_names / i) {
             if (k >= 0) {
-                res[node_names[i]] =  node_names[k]; 
+                if (node_names / k) {
+                    res[node_names[i]] =  node_names[k]; 
+                } else {
+                    res[node_names[i]] =  "root";
+                }
             } else {
                 res[node_names[i]] = None;
             }

src/core/global_things.cpp

@@ -115,7 +115,7 @@ _String const kEmptyString, kPromptForFilePlaceholder("PROMPT_FOR_FILE"),
                     "\"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.85"),
+    kHyPhyVersion = _String("2.5.86"),
 
     kNoneToken = "None", kNullToken = "null",
     kNoKWMatch = "__input_value_not_given__",

src/core/matrix.cpp

@@ -4474,7 +4474,7 @@ void _Matrix::Multiply(_Matrix &storage, _Matrix const &secondArg) const
                 __m256d value_op = _mm256_set1_pd(value);
 
 #ifdef _SLKP_USE_FMA3_INTRINSICS
-#define CELL_OP(x)                                                             \
+#define CELL_OP_2(x)                                                           \
   _mm256_storeu_pd(res + x,                                                    \
                    _mm256_fmadd_pd(value_op, _mm256_loadu_pd(secArg + x),      \
                                    _mm256_loadu_pd(res + x)))

src/core/topology.cpp

@@ -747,6 +747,16 @@ _List _TreeTopology::RemoveANode(HBLObjectRef nodeName) {
       } else {
         throw _String("Matrix-valued argument was expected to contain strings");
       }
+    } else if (nodeName->ObjectClass() == ASSOCIATIVE_LIST) {
+      _AssociativeList *names = (_AssociativeList *)nodeName;
+      for (AVLListXLIteratorKeyValue key_value : names->ListIterator()) {
+        HBLObjectRef var_value = (HBLObjectRef)key_value.get_object();
+        if (var_value->ObjectClass() == STRING) {
+          RemoveNodeByName((_FString *)var_value);
+        }
+      }
+      clean_indices.Sort();
+      _RemoveNodeList(clean_indices);
     } else {
       throw _String(
           "An invalid argument (not a string or a string matrix) supplied");

src/mains/unix.cpp

@@ -254,10 +254,17 @@ _String getLibraryPath() {
   _String baseDir(buffer);
   baseDir.Trim(0, baseDir.FindBackwards(dirSlash, 0, -1) - 1L);
 #else
-  char curWd[_HYPHY_MAX_PATH_LENGTH];
-  getcwd(curWd, _HYPHY_MAX_PATH_LENGTH);
+  char *curWd = (char *)malloc(_HYPHY_MAX_PATH_LENGTH);
+  if (!curWd || !getcwd(curWd, _HYPHY_MAX_PATH_LENGTH)) {
+    fprintf(stderr,
+            "Failed to get current working directory path; possibly because it "
+            "is longer than %ld characters",
+            _HYPHY_MAX_PATH_LENGTH);
+    exit(1);
+  }
 
   _String baseDir(curWd);
+  free(curWd);
 #endif
 
   ensure_trailing_dirSlash(baseDir);
@@ -814,12 +821,20 @@ int main(int argc, char *argv[]) {
   printf ("%ld\n", sscanf (" 0.1e2 beavis", "%lf%n", &value, &read));
   */
 
-  char curWd[4096], dirSlash = get_platform_directory_char();
+  char *curWd = (char *)malloc(4096), dirSlash = get_platform_directory_char();
 
-  getcwd(curWd, 4096);
+  if (!curWd || !getcwd(curWd, 4096)) { // failed to get cwd
+    fprintf(stderr,
+            "Failed to get current working directory path; possibly because it "
+            "is longer than %d characters",
+            4096);
+    return 1;
+  }
 
   _String baseDir(curWd);
 
+  free(curWd);
+
   auto ensure_trailing_dirSlash = [&dirSlash](auto &str) {
     if (str.get_char(str.length() - 1) != dirSlash) {
       str = str & dirSlash;