Merge branch 'develop'

CMakeLists.txt

@@ -495,15 +495,7 @@ else(${OPENMP_FOUND})
     target_link_libraries(HYPHYDEBUG ${DEFAULT_LIBRARIES})
 endif(${OPENMP_FOUND})
 
-if(${MPI_FOUND})
-    set_property(
-        TARGET HYPHYDEBUG
-        APPEND PROPERTY COMPILE_DEFINITIONS __HYPHYMPI__
-    )
-    target_link_libraries(HYPHYDEBUG ${DEFAULT_LIBRARIES} ${MPI_LIBRARIES} OpenMP::OpenMP_CXX)
-else(${MPI_FOUND})
     target_link_libraries(HYPHYDEBUG ${DEFAULT_LIBRARIES})
-endif(${MPI_FOUND})
 
 
 add_custom_target(DEBUG DEPENDS HYPHYDEBUG)

res/TemplateBatchFiles/GARD.bf

@@ -201,7 +201,10 @@ utility.ForEach (alignments.Extract_site_patterns ("gard.filter"), "_pattern_",
 
 gard.variableSiteMap = Transpose (utility.DictToArray (gard.variableSiteMap)) % 0; // sort by 1st column
 gard.variableSites = Rows (gard.variableSiteMap);
-gard.inverseVariableSiteMap = unility.SwapKeysAndValues(gard.variableSiteMap);
+gard.inverseVariableSiteMap = {};
+for (index, pattern; in; gard.variableSiteMap) {
+    gard.inverseVariableSiteMap[pattern] = index;
+}
 gard.numberOfPotentialBreakPoints = gard.variableSites - 1;
 
 io.ReportProgressMessage ("", ">Loaded a `gard.dataType` multiple sequence alignment with **`gard.numSeqs`** sequences, **`gard.numSites`** sites (`gard.variableSites` of which are variable) from \`" +
@@ -606,21 +609,17 @@ lfunction gard.createLikelihoodFunctionForExport (id,  model) {
 */
 lfunction gard.validatePartititon (definition, minSize, totalSites) {
     lastBP = 0;
-    //console.log ("");
-    //console.log (definition);
-    bpCount = utility.Array1D (definition);
-    for (i = 0; i < bpCount; i+=1) {
-        if (definition[i] - lastBP + 1 < minSize) {
-            //console.log ("FAIL");
+
+    for (current_bp; in; definition) {
+        if (current_bp - lastBP + 1 < minSize) {
             return FALSE;
         }
-        lastBP = definition[i];
+        lastBP = current_bp;        
     }
+
     if (totalSites - lastBP < minSize) {
-        //console.log ("FAIL");
         return FALSE;
     }
-    //console.log ("PASS");
     return TRUE;
 }
 
@@ -631,7 +630,8 @@ lfunction gard.validatePartititon (definition, minSize, totalSites) {
 
  */
 lfunction gard.modelIsNotInMasterList(masterList, breakPoints) {
-    return utility.KeyExists(masterList, '' + breakPoints) == FALSE;
+    return masterList / ('' + breakPoints) == FALSE;
+    //return utility.KeyExists(masterList, '' + breakPoints) == FALSE;
 }
 
 function gard.concludeAnalysis(bestOverallModel) {
@@ -807,7 +807,7 @@ lfunction gard.GA.recombineModels (parentModels, populationSize) {
 
             breakPoints = {numberOfBreakPoints, 1};
             for (breakPointNumber=0; breakPointNumber < numberOfBreakPoints; breakPointNumber += 1) {
-                if (random.TRUE_or_FALSE()) {
+                if (Random(0,1) <= 0.5) {
                     breakPoints[breakPointNumber] = parentModel1[breakPointNumber];
                 } else {
                     breakPoints[breakPointNumber] = parentModel2[breakPointNumber];
@@ -835,16 +835,10 @@ lfunction gard.GA.recombineModels (parentModels, populationSize) {
  * @returns nothing... the gard.GA.storeMultiBreakPointModelResults function gets called which updates the intergenerationalModels object
  */
 function gard.GA.evaluateModels (models) {
-    modelIds = utility.Keys(models);
-    numberOfModels = Columns(modelIds);
-
-    for(modelIndex=0; modelIndex<numberOfModels; modelIndex += 1) {
-        modelId = modelIds[modelIndex];
-        cAIC = models[modelId];
-
-
+
+    for (modelID, cAIC; in; models)  {
         if (cAIC == math.Infinity) {
-            breakPoints = gard.Helper.convertMatrixStringToMatrix(modelId);
+            breakPoints = gard.Helper.convertMatrixStringToMatrix(modelID);
             mpi.QueueJob (gard.queue, "gard.obtainModel_cAIC", {"0" : breakPoints__,
                                                      "1" : gard.model,
                                                      "2" : gard.baseLikelihoodInfo},
@@ -900,8 +894,7 @@ lfunction gard.GA.selectModels (evaluatedModels, numberOfModelsToKeep) {
  * @returns a {Bolean}
  */
 lfunction gard.GA.modelSetsAreTheSame(modelSet1, modelSet2) {
-    s1 = modelSet1;
-    return (s1-modelSet2) == 0;
+    return modelSet1 == modelSet2;
 }
 
 /**
@@ -930,20 +923,19 @@ lfunction gard.GA.generateNewGenerationOfModelsByMutatingModelSet(parentModels,
         while(modelIsValid == FALSE && failedAttempts < ^"gard.maxFailedAttemptsToMakeNewModel") {
             parentModel = gard.Helper.convertMatrixStringToMatrix(modelIds[i]);
             breakPoints = {numberOfBreakPoints, 1};
-            for(breakPointIndex=0; breakPointIndex<numberOfBreakPoints; breakPointIndex=breakPointIndex+1) {
+            for(breakPointIndex=0; breakPointIndex<numberOfBreakPoints; breakPointIndex += 1) {
 
                 if(Random(0,1) < mutationRate) { // keep the break point the same
                     breakPoints[breakPointIndex] = parentModel[breakPointIndex];
                 } else {
                     if(Random(0,1) < rateOfMutationsThatAreSmallShifts) { // move the break point by a random small amount
                         notValid = TRUE;
                         while (notValid) {
-                            distanceOfStep = Min (1,random.poisson(2));
-                            if (random.TRUE_or_FALSE()) { // randomly decide if the break point moves right or left
-                                distanceOfStep = - distanceOfStep;
+                            distanceOfStep = Max (1,random.poisson(1));
+                            if (Random (0,1) <= 0.5) { // randomly decide if the break point moves right or left
+                                distanceOfStep = -distanceOfStep;
                             }
-                            variableSiteMapIndexOfParentBreakPoint = utility.Find(^"gard.variableSiteMap", parentModel[breakPointIndex]);
-                            variableSiteMapIndexOfParentBreakPoint += distanceOfStep;
+                            variableSiteMapIndexOfParentBreakPoint = (^"gard.inverseVariableSiteMap")[parentModel[breakPointIndex]] + distanceOfStep;
                             if (variableSiteMapIndexOfParentBreakPoint >= 0 && variableSiteMapIndexOfParentBreakPoint < (^"gard.variableSites")) {
                                 newBreakPoint = (^"gard.variableSiteMap")[variableSiteMapIndexOfParentBreakPoint];
                                 notValid = FALSE;
@@ -1019,14 +1011,4 @@ lfunction gard.Helper.convertMatrixStringToMatrix(matrixString){
 //sort numeric Nx1 matrix into assending order.
 lfunction gard.Helper.sortedMatrix(matrix) {
     return matrix % 0;
-
-    /*sortedMatrix = {1,Columns(matrix)};
-
-    for (i=0; i<Columns(matrix); i+=1) {
-        minElementLeft=Min(matrix, 1);
-        sortedMatrix[i] = minElementLeft[0];
-        matrix[minElementLeft[1]] = ^"math.Infinity";
-    }
-
-    return sortedMatrix;*/
 }

res/TemplateBatchFiles/SelectionAnalyses/MEME.bf

@@ -508,7 +508,7 @@ lfunction meme.handle_a_site (lf_fel, lf_bsrel, filter_data, partition_index, pa
     //io.SpoolLF (lf_bsrel, "/tmp/meme.debug", "MEME");
 
     Optimize (results, ^lf_bsrel, {
-            //"OPTIMIZATION_METHOD" : "nedler-mead",
+            //"OPTIMIZATION_METHOD" : "gradient-descent",
             "OPTIMIZATION_START_GRID" : 
              {
                 "0" : {

res/TemplateBatchFiles/files.lst

@@ -10,6 +10,7 @@
 
 "","A collection of tools for evolutionary hypothesis testing coding-sequence data.","!Evolutionary Hypothesis Testing";
 "Contrast-FEL","Perform a site-level test for differences in selective pressures between predetermined sets of branches.","SelectionAnalyses/FEL-contrast.bf";
+"FMM", "Fit a model that permits double (and triple) instantaneous nucleotide substitutions","SelectionAnalyses/FitMultiModel.bf";
 
 "","Evolutionary rates on non-coding data.","!Relative evolutionary rate inference";
 "LEISR","Infer relative evolutionary rates on a nucleotide or protein alignment, in a spirit similar to Rate4Site (PMID: 12169533).","LEISR.bf";

res/TemplateBatchFiles/libv3/UtilityFunctions.bf

@@ -142,9 +142,8 @@ lfunction utility.CheckCacheFile (data_info) {
  * @returns -1 if value is not found, otherwise returns the position
  */
 lfunction utility.Find (array, value) {
-    d = Rows (array) * Columns (array);
-    for (i = 0; i < d; i+=1) {
-        if (array [i] == value) {
+    for (i, v; in; array) {
+        if (v == value) {
             return i;
         }
     }

res/TemplateBatchFiles/libv3/all-terms.bf

@@ -116,6 +116,11 @@ namespace terms{
     function aminoacidRate(fromA, toA) {
         return "Substitution rate from amino-acid " + fromA + " to amino-acid " + toA;
     }
+
+    function propertyImportance  (id) {
+        return "Importance factor for property " + id;
+    }
+
     function binaryRate(fromX, toX) {
         return "Substitution rate from character " + fromX + " to character " + toX;
     }
@@ -368,6 +373,7 @@ namespace terms{
         length_expression       = "length-expression";
 
         data                    = "data";
+        residue_properties      = "residue_properties";
 
     }
 
@@ -401,6 +407,7 @@ namespace terms{
         synonymous_rate               = "synonymous rate";
         nonsynonymous_rate            = "non-synonymous rate";
         omega_ratio                   = "non-synonymous/synonymous rate ratio";
+        log_omega_ratio               = "log of non-synonymous/synonymous rate ratio";
         multiple_hit_rate             = "rate at which 2 nucleotides are changed instantly within a single codon";
         triple_hit_rate               = "rate at which 3 nucleotides are changed instantly within a single codon";
         triple_hit_rate_syn           = "rate at which 3 nucleotides are changed instantly within a single codon between synonymous codon islands";

res/TemplateBatchFiles/libv3/models/parameters.bf

@@ -824,8 +824,31 @@ lfunction parameters.SetCategoryVariables (model) {
     }
 
     return cat_var_substitutions;
-
-
-
+}
+
+/**
+ * Given a set of strings, convert them to valid IDs that don't clash following renames
+ * @name parameters.ValidateID
+ * @param {Array/Dict} ids - the ids to validate (should be unique), if dict, should be i -> ID
+ * @returns {Dictionary} a dictionary containing old -> new id map
+*/
+
+lfunction parameters.ValidateIDs (ids) {
+    result = {};
+    N = utility.Array1D (ids);
+    for (i = 0; i < N; i+=1) {
+        try_name = ids [i] && 7;
+        if (result / try_name ) {
+            for (k = 2; ; k+=1) {
+                tagged = try_name + "_" + k;
+                if (!(result/tagged)) {
+                    try_name = tagged;
+                    break;
+                }
+            }
+        }
+        result [ids[i]] = try_name;
+    }
+    return result;
 }
 

res/TemplateBatchFiles/libv3/tasks/distances.bf

@@ -23,7 +23,7 @@ lfunction distances.nucleotide.tn93 (filter, freqs, options) {
 		K2P = TRUE;
 	}
 	else {
-	    K2P = FASLE;
+	    K2P = FALSE;
 		K1 = 2*freqs[0]*freqs[2]/fR;
 		K2 = 2*freqs[1]*freqs[3]/fY;
 		K3 = 2*(fR*fY-freqs[0]*freqs[2]*fY/fR-freqs[1]*freqs[3]*fR/fY);

res/TemplateBatchFiles/libv3/tasks/mpi.bf

@@ -180,6 +180,8 @@ namespace mpi {
 
 
                     send_to_nodes * 0;
+                    queue ["cache"] = {};
+
                 }
             }
 
@@ -197,6 +199,8 @@ namespace mpi {
             for (k = 1; k < mpi_node_count; k += 1) {
                 queue [k] = {utility.getGlobalValue("terms.mpi.job_id") : None, utility.getGlobalValue("terms.mpi.callback") : None, utility.getGlobalValue("terms.mpi.arguments"): None};
             }
+
+            // this will store jobs that were previously sent to each node; avoiding redefinition if possible
         }
         return queue;
     }
@@ -225,8 +229,21 @@ namespace mpi {
                 node = aux._handle_receieve (queue);
             }
 
-
-            complete_function_dump = aux.queue_export_function (job);
+            complete_function_dump = None;
+
+            if (utility.Has (queue, "cache" , "AssociativeList")) {
+                if ((queue["cache"])[node] == job) {
+                    complete_function_dump = "";
+                    //console.log ("CACHED MPI preamble"); 
+                } else {
+                    (queue["cache"])[node] = job;
+                }
+            } 
+
+            if (None == complete_function_dump) {
+                complete_function_dump = aux.queue_export_function (job);
+            }
+
             //console.log (complete_function_dump);
             job_id = get_next_job_id();
             //fprintf (stdout, "Sending to node ", node, "\n");
@@ -261,6 +278,8 @@ namespace mpi {
                 }
             } while (node < mpi_node_count);
         }
+
+        queue = None
     }
 
     namespace aux {

res/TemplateBatchFiles/libv3/tasks/trees.bf

@@ -428,7 +428,10 @@ lfunction trees.ExtractTreeInfo(tree_string) {
 
     rooted = utility.Array1D ((flat_tree[(flat_tree[0])["Root"]])["Children"]) == 2;
 
-    DeleteObject (flat_tree, branch_lengths, branch_names, branch_count);
+    flat_tree       = None;
+    branch_lengths  = None;
+    branch_names    = None;
+    branch_count    = None;
 
     return {
         ^"terms.trees.newick": Format(T, 1, 0),
@@ -1093,12 +1096,12 @@ lfunction tree.infer.NJ (datafilter, distances) {
                 treeNodes[i][2] = njm[i][2]; // branch length
             }
 
-            DeleteObject (njm);
+            njm = None;
 		}
 	}
 
     if (flush_distances) {
-        DeleteObject (distances);
+        distances = None;
     }
 
     return tree._matrix2string (treeNodes, N, alignments.GetSequenceNames (datafilter), TRUE);

src/core/avllistxl_iterator.cpp

@@ -39,8 +39,6 @@
 
 #include "avllistxl_iterator.h"
 
-#define  AVL_LISTXL_ITERATOR_ENDINDEX (-1L)
-
 AVLListXLIterator::AVLListXLIterator (_AVLListXL const * the_list):
   current_location (AVL_LISTXL_ITERATOR_ENDINDEX),
   traversal_lookup (AVL_LISTXL_ITERATOR_ENDINDEX),