diff --git a/res/TemplateBatchFiles/SelectionAnalyses/RELAX.bf b/res/TemplateBatchFiles/SelectionAnalyses/RELAX.bf
index 7597a31e2..64b64c652 100644
--- a/res/TemplateBatchFiles/SelectionAnalyses/RELAX.bf
+++ b/res/TemplateBatchFiles/SelectionAnalyses/RELAX.bf
@@ -283,21 +283,35 @@ if (relax.model_set == "All") { // run all the models
                 relax.filter_names,
                 None);
 
+        relax.distribution          = models.codon.BS_REL.ExtractMixtureDistribution(relax.ge.bsrel_model);
+        relax.weight_multipliers    = parameters.helper.stick_breaking (utility.SwapKeysAndValues(utility.MatrixToDict(relax.distribution["weights"])),None);
+        relax.constrain_parameters   = parameters.ConstrainMeanOfSet(relax.distribution["rates"],relax.weight_multipliers,1,"relax");
+        
+        for (key, value; in; relax.constrain_parameters[terms.global]){
+            model.generic.AddGlobal (relax.ge.bsrel_model, value, key);
+            parameters.SetRange (value, terms.range_almost_01);
+        }
+        
+        relax.distribution["rates"] = Transpose (utility.Values (relax.constrain_parameters[terms.global]));
+        
         for (relax.i = 1; relax.i < relax.rate_classes; relax.i += 1) {
             parameters.SetRange (model.generic.GetGlobalParameter (relax.ge.bsrel_model , terms.AddCategory (terms.parameters.omega_ratio,relax.i)), terms.range_almost_01);
         }
         parameters.SetRange (model.generic.GetGlobalParameter (relax.ge.bsrel_model , terms.AddCategory (terms.parameters.omega_ratio,relax.rate_classes)), terms.range_gte1);
         
+        // constrain the mean of this distribution to 1
+        
 
         relax.model_object_map = { "relax.ge" :       relax.ge.bsrel_model };
 
         io.ReportProgressMessageMD ("RELAX", "gd", "Fitting the general descriptive (separate k per branch) model");
         selection.io.startTimer (relax.json [terms.json.timers], "General descriptive model fitting", 2);
 
-        relax.distribution = models.codon.BS_REL.ExtractMixtureDistribution(relax.ge.bsrel_model);
         PARAMETER_GROUPING = {};
         PARAMETER_GROUPING + relax.distribution["rates"];
         PARAMETER_GROUPING + relax.distribution["weights"];
+        
+        
 
         if (Type (relax.ge_guess) != "Matrix") {
             // first time in 
@@ -315,8 +329,6 @@ if (relax.model_set == "All") { // run all the models
             
             
             parameters.DeclareGlobalWithRanges ("relax.bl.scaler", 1, 0, 1000);
-                        
-            //VERBOSITY_LEVEL = 10;             
                          
             relax.grid_search.results =  estimators.FitLF (relax.filter_names, relax.trees,{ "0" : {"DEFAULT" : "relax.ge"}},
                                         relax.final_partitioned_mg_results,
@@ -430,6 +442,12 @@ if (relax.model_set == "All") { // run all the models
                                                      0,
                                                      relax.k_estimates);
 
+
+        for (relax.i = 1; relax.i <= relax.rate_classes; relax.i += 1) {
+            //console.log (model.generic.GetGlobalParameter (relax.ge.bsrel_model , terms.AddCategory (terms.parameters.omega_ratio,relax.i)));
+            parameters.RemoveConstraint (model.generic.GetGlobalParameter (relax.ge.bsrel_model , terms.AddCategory (terms.parameters.omega_ratio,relax.i)));
+        }
+
         break;
     }
 
diff --git a/res/TemplateBatchFiles/SelectionAnalyses/SLAC.bf b/res/TemplateBatchFiles/SelectionAnalyses/SLAC.bf
index c10e0d0e4..2ed782984 100644
--- a/res/TemplateBatchFiles/SelectionAnalyses/SLAC.bf
+++ b/res/TemplateBatchFiles/SelectionAnalyses/SLAC.bf
@@ -517,6 +517,7 @@ lfunction slac.compute_the_counts (matrix, tree, lookup, selected_branches, coun
 
 
     selected_branch_total_length  = +selected_branches_lengths;
+        
     io.CheckAssertion ("`&selected_branch_total_length`>0", "SLAC cannot be applied to a branch selection with total zero branch length (i.e. no variation)");
 
     /*  columns
@@ -593,7 +594,7 @@ lfunction slac.compute_the_counts (matrix, tree, lookup, selected_branches, coun
                     relative_branch_length = 1;
                     fr = Eval (fully_resolved);
                     relative_branch_length = selected_branches_lengths[i] / selected_branch_total_length;
-
+                    
                     for (k = 0; k < 4; k += 1) {
                         report_resolved_by_branch[i*column_count + k] += fr[k];
                         report_averaged_by_branch[i*column_count + k] += fr[k];
@@ -605,12 +606,12 @@ lfunction slac.compute_the_counts (matrix, tree, lookup, selected_branches, coun
                             by_site_scaler [s] += (-selected_branches_lengths[i]);
                             psi = parent_state*state_count+parent_state;
 
-
-                            fr = Eval (fully_resolved);
+                            /*fr = Eval (fully_resolved);
+                            
                             for (k = 0; k < 2; k += 1) {
                                 report_averaged[s*column_count + k] += fr[k];
                                 report_resolved[s*column_count + k] += fr[k];
-                            }
+                            }*/
 
                             relative_branch_length = 1;
                             fr = Eval (fully_resolved);
diff --git a/res/TemplateBatchFiles/SelectionAnalyses/contrast-fel.bf b/res/TemplateBatchFiles/SelectionAnalyses/contrast-fel.bf
index f39d096ef..0e7737472 100644
--- a/res/TemplateBatchFiles/SelectionAnalyses/contrast-fel.bf
+++ b/res/TemplateBatchFiles/SelectionAnalyses/contrast-fel.bf
@@ -98,7 +98,8 @@ selection.io.startTimer (fel.json [terms.json.timers], "Total time", 0);
 
 namespace fel {
     LoadFunctionLibrary ("modules/shared-load-file.bf");
-    load_file ({utility.getGlobalValue("terms.prefix"): "fel", utility.getGlobalValue("terms.settings") : {utility.getGlobalValue("terms.settings.branch_selector") : "fel.select_branches"}});
+    load_file ({utility.getGlobalValue("terms.prefix"): "fel", 
+                utility.getGlobalValue("terms.settings") : {utility.getGlobalValue("terms.settings.branch_selector") : "fel.select_branches"}});
 }
 
 KeywordArgument ("srv", "Include synonymous rate variation in the model", "Yes");
diff --git a/res/TemplateBatchFiles/SelectionAnalyses/modules/io_functions.ibf b/res/TemplateBatchFiles/SelectionAnalyses/modules/io_functions.ibf
index bd4c1bf0d..5f810711d 100644
--- a/res/TemplateBatchFiles/SelectionAnalyses/modules/io_functions.ibf
+++ b/res/TemplateBatchFiles/SelectionAnalyses/modules/io_functions.ibf
@@ -19,14 +19,14 @@ lfunction selection.io.defineBranchSets(partition_info) {
         tree_for_analysis = (partition_info[k])[utility.getGlobalValue("terms.data.tree")];
         available_models = {};
         utility.ForEach (tree_for_analysis[utility.getGlobalValue("terms.trees.model_map")], "_value_", "`&available_models`[_value_] += 1");
-        
+
         for (k2; in; utility.Keys   (available_models)) {
             key_counts[k2] = key_counts[k2] + 1;
         }
     }
-    
-    
-    
+
+
+
     list_models = {};
     for (k,v; in; key_counts) {
         if (v == Abs (partition_info)) {
@@ -39,11 +39,11 @@ lfunction selection.io.defineBranchSets(partition_info) {
         list_models   = utility.sortStrings(utility.Keys   (list_models));
     }
 
-     
+
     selectTheseForTesting = {
         option_count + 3, 2
     };
-    
+
 
     selectTheseForTesting[0][0] = "All";
     selectTheseForTesting[0][1] = "Include all branches in the analysis";
@@ -61,7 +61,7 @@ lfunction selection.io.defineBranchSets(partition_info) {
             selectTheseForTesting[3 + k][1] = "Set of " + available_models[list_models[k]] + " unlabeled branches";
         }
     }
-    
+
     ChoiceList(testSet, "Choose the set of branches to test for selection", 1, NO_SKIP, selectTheseForTesting);
 
     io.CheckAssertion ("`&testSet` >= 0", "User cancelled branch selection; analysis terminating");
@@ -195,17 +195,17 @@ lfunction selection.io.report_viterbi_path (path) {
     rates = Rows (distribution);
     settings ["header"] = FALSE;
     last_state = path[0];
-    
-    
+
+
     row_matrix = {{index, last_state, rate}};
-    
+
     for (index, rate; in; path) {
         if (rate != last_state) {
             fprintf (stdout, io.FormatTableRow (row_matrix, settings));
             last_state = rate;
         }
     }
-    
+
     fprintf (stdout, "\n");
     return row_matrix;
 }
@@ -414,6 +414,15 @@ function selection.io.json_store_branch_attribute(json, attribute_name, attribut
     ((json[terms.json.branch_attributes])[terms.json.attribute])[attribute_name] = {terms.json.attribute_type : attribute_type,
                                                                                             terms.json.display_order: display_order};
 
+
+
+    for (selection.io.json_store_branch_attribute.branch_name,selection.io.json_store_branch_attribute.branch_tag; in; values) {
+      utility.EnsureKey ((json[terms.json.branch_attributes])[partition], selection.io.json_store_branch_attribute.branch_name);
+      (((json[terms.json.branch_attributes])[partition])[selection.io.json_store_branch_attribute.branch_name])[attribute_name] = selection.io.json_store_branch_attribute.branch_tag;
+    }
+
+
+    /*
     utility.ForEach (utility.Keys (values), "selection.io.json_store_branch_attribute.branch_name",
                              "utility.EnsureKey ((json[terms.json.branch_attributes])[partition], selection.io.json_store_branch_attribute.branch_name)");
 
@@ -422,6 +431,7 @@ function selection.io.json_store_branch_attribute(json, attribute_name, attribut
                                     (((json[terms.json.branch_attributes])[partition])[selection.io.json_store_branch_attribute.branch_name])[attribute_name] = values[selection.io.json_store_branch_attribute.branch_name];
                                 }
                               ");
+    */
  }
 
 
@@ -436,4 +446,3 @@ function selection.io.extract_branch_info(branch_spec, callback) {
     branch_spec["selection.io._aux.extract_branch_info.callback"][""];
     return selection.io.extract_branch_info_result;
 }
-
diff --git a/res/TemplateBatchFiles/SelectionAnalyses/modules/shared-load-file.bf b/res/TemplateBatchFiles/SelectionAnalyses/modules/shared-load-file.bf
index 017417ae8..b6c12f654 100644
--- a/res/TemplateBatchFiles/SelectionAnalyses/modules/shared-load-file.bf
+++ b/res/TemplateBatchFiles/SelectionAnalyses/modules/shared-load-file.bf
@@ -72,6 +72,7 @@ function load_file (prefix) {
 
     sample_size=codon_data_info[utility.getGlobalValue("terms.data.sites")]*codon_data_info[utility.getGlobalValue("terms.data.sequences")];
 
+
     codon_data_info[utility.getGlobalValue("terms.data.sample_size")] = sample_size;
     upper_prefix = prefix && 1; //uppercase the prefix for json name
     codon_data_info[utility.getGlobalValue("terms.json.json")] = codon_data_info[utility.getGlobalValue("terms.data.file")] + "."+upper_prefix+".json";
@@ -86,8 +87,19 @@ function load_file (prefix) {
         name_mapping = {};
         utility.ForEach (alignments.GetSequenceNames (prefix+".codon_data"), "_value_", "`&name_mapping`[_value_] = _value_");
     }
-    
-    
+
+    // check for duplicates
+    duplicate_sequences = codon_data_info[^"terms.data.sequences"] - alignments.HasDuplicateSequences (codon_data_info[^"terms.data.datafilter"],-1);
+    if (duplicate_sequences > 0) {
+      fprintf(stdout, "\n-------\n", io.FormatLongStringToWidth(
+      ">[WARNING] This dataset contains " + duplicate_sequences + " duplicate " + io.SingularOrPlural (duplicate_sequences, 'sequence', 'sequences') + ".
+      Identical sequences do not contribute any information to the analysis and only slow down computation.
+      Please consider removing duplicate or 'nearly' duplicate sequences,
+      e.g. using https://github.com/veg/hyphy-analyses/tree/master/remove-duplicates
+      prior to running selection analyses", 72),
+      "\n-------\n");
+    }
+
     utility.SetEnvVariable(utility.getGlobalValue ("terms.trees.data_for_neighbor_joining"),
                            codon_data_info[utility.getGlobalValue("terms.data.datafilter")]);
 
@@ -122,7 +134,7 @@ function load_file (prefix) {
             },
             ...
         */
-        
+
 
     partition_count = Abs (partitions_and_trees);
 
@@ -203,7 +215,7 @@ function load_file (prefix) {
 
 function store_tree_information () {
     // Place in own attribute called `tested`
-    
+
      selection.io.json_store_key_value_pair (json, None, utility.getGlobalValue("terms.json.tested"), selected_branches);
 
         /**  this will return a dictionary of selected branches; one set per partition, like in
@@ -277,30 +289,30 @@ function doGTR (prefix) {
                                         terms.nucleotideRate ("G","T") : { utility.getGlobalValue ("terms.fit.MLE") : 0.25}
                                     }
                                  });
-            
+
 
     //utility.ToggleEnvVariable("VERBOSITY_LEVEL", 10);
 
     gtr_results = estimators.FitGTR(filter_names,
                                          trees,
                                          gtr_results);
-                                         
-                                        
+
+
     KeywordArgument ("kill-zero-lengths", "Automatically delete internal zero-length branches for computational efficiency (will not affect results otherwise)", "Yes");
-                       
+
     kill0 = io.SelectAnOption (
         {
             "Yes":"Automatically delete internal zero-length branches for computational efficiency (will not affect results otherwise)",
             "No":"Keep all branches"
-        }, 
+        },
         "The set of properties to use in the model") == "Yes";
-        
-         
-    if (kill0) {                   
+
+
+    if (kill0) {
         for (index, tree; in; trees) {
             deleted = {};
             if (^(prefix + ".selected_branches") / index) {
-                trees[index] = trees.KillZeroBranches (tree, (gtr_results[^"terms.branch_length"])[index], (^(prefix + ".selected_branches"))[index], deleted);       
+                trees[index] = trees.KillZeroBranches (tree, (gtr_results[^"terms.branch_length"])[index], (^(prefix + ".selected_branches"))[index], deleted);
             } else {
                 trees[index] = trees.KillZeroBranches (tree, (gtr_results[^"terms.branch_length"])[index], null, deleted);
             }
@@ -314,11 +326,11 @@ function doGTR (prefix) {
         }
         for (i = 0; i < partition_count; i+=1) {
             (partitions_and_trees[i])[^"terms.data.tree"] = trees[i];
-        }   
+        }
         store_tree_information ();
     }
-    
-                                         
+
+
     io.ReportProgressMessageMD (prefix, "nuc-fit", "* " +
         selection.io.report_fit (gtr_results, 0, 3*(^"`prefix`.sample_size")));
 
@@ -393,7 +405,7 @@ function doPartitionedMG (prefix, keep_lf) {
 
     if (partition_count > 1) {
        //partition_scalers = selection.io.extract_global_MLE_re (partitioned_mg_results, "^" + utility.getGlobalValue("terms.parameters.omega_ratio"));
-   
+
     }
 
     /** extract and report dN/dS estimates */
diff --git a/res/TemplateBatchFiles/TemplateModels/chooseGeneticCode.def b/res/TemplateBatchFiles/TemplateModels/chooseGeneticCode.def
index f106757f8..a0a4877c4 100644
--- a/res/TemplateBatchFiles/TemplateModels/chooseGeneticCode.def
+++ b/res/TemplateBatchFiles/TemplateModels/chooseGeneticCode.def
@@ -540,21 +540,14 @@ lfunction CompareCodonProperties(codon1, codon2, code)
 /*----------------------------------------------------------------------------------------------------------*/
 
 function defineCodonToAA() {
-    codonToAAMap = {};
-    nucChars = "ACGT";
-
-    for (p1 = 0; p1 < 64; p1 = p1 + 1) {
-        codonToAAMap[nucChars[p1$16] + nucChars[p1 % 16 $4] + nucChars[p1 % 4]] = _hyphyAAOrdering[_Genetic_Code[p1]];
-    }
-
-    return codonToAAMap;
+    return defineCodonToAAGivenCode(_Genetic_Code);
 }
 
 /*----------------------------------------------------------------------------------------------------------*/
 
 function defineCodonToAAGivenCode(code) {
     codonToAAMap = {};
-    nucChars = "ACGT";
+    nucChars     = "ACGT";
 
     for (p1 = 0; p1 < 64; p1 += 1) {
         codonToAAMap[nucChars[p1$16] + nucChars[p1 % 16 $4] + nucChars[p1 % 4]] = _hyphyAAOrdering[code[p1]];
diff --git a/res/TemplateBatchFiles/files.lst b/res/TemplateBatchFiles/files.lst
index 28be624cc..534d1cee7 100644
--- a/res/TemplateBatchFiles/files.lst
+++ b/res/TemplateBatchFiles/files.lst
@@ -95,7 +95,6 @@
 
 "","Test for recombination.","!Recombination";
 "GARD","[GARD] Screen an alignment using GARD (requires an MPI environment).","GARD.bf";
-"GRDR","Process GARD results.","GARDProcessor.bf";
 "LHT","A Likelihood Ratio Test to detect conflicting phylogenetic signal Huelsenbeck and Bull, 1996. [Contributed by Olivier Fedrigo].","LHT.bf";
 "SBL","Search an alignment for a single breakpoint.","SingleBreakpointRecomb.bf";
 "SPL","Plot genetic distances (similarity) of one sequence against all others in an alignment, using a sliding window. Optionally, determine NJ-based clustering and bootstrap support in every window. This is a HyPhy adaptation of the excellent (but Windows only tool) SimPlot (and/or VarPlot) written by Stuart Ray (http://sray.med.som.jhmi.edu/SCRoftware/simplot/)","SimilarityPlot.bf";
diff --git a/res/TemplateBatchFiles/libv3/UtilityFunctions.bf b/res/TemplateBatchFiles/libv3/UtilityFunctions.bf
index a7ef78a67..587459c6b 100644
--- a/res/TemplateBatchFiles/libv3/UtilityFunctions.bf
+++ b/res/TemplateBatchFiles/libv3/UtilityFunctions.bf
@@ -298,7 +298,7 @@ function utility.MapWithKey (object, key_name, lambda_name, transform) {
             utility.MapWithKey.return_object = {utility.MapWithKey.rows,  utility.MapWithKey.columns};
 
             ^(lambda_name) := object [utility.MapWithKey.r][utility.MapWithKey.c];
-            ^(key_name) := {{utility.MapWithKey.r,utility.MapWithKey.c}}
+            ^(key_name) := {{utility.MapWithKey.r,utility.MapWithKey.c}};
             for (utility.MapWithKey.r = 0; utility.MapWithKey.r < utility.MapWithKey.rows; utility.MapWithKey.r += 1) {
                 for (utility.MapWithKey.c = 0; utility.MapWithKey.c < utility.MapWithKey.columns; utility.MapWithKey.c += 1) {
                     utility.MapWithKey.temp = Eval (transform);
diff --git a/res/TemplateBatchFiles/libv3/models/codon/MG_REV.bf b/res/TemplateBatchFiles/libv3/models/codon/MG_REV.bf
index 6b801c6cd..f89ceb3ce 100644
--- a/res/TemplateBatchFiles/libv3/models/codon/MG_REV.bf
+++ b/res/TemplateBatchFiles/libv3/models/codon/MG_REV.bf
@@ -14,7 +14,7 @@ LoadFunctionLibrary("../../UtilityFunctions.bf");
 lfunction models.codon.MG_REV.ModelDescription(type, code) {
 
     codons = models.codon.MapCode(code);
-    
+
     return {
         utility.getGlobalValue("terms.alphabet"): codons[utility.getGlobalValue("terms.sense_codons")],
         utility.getGlobalValue("terms.bases"): utility.getGlobalValue("models.DNA.alphabet"),
@@ -122,8 +122,8 @@ 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
-          
-        if (Type(value) == "AssociativeList") { 
+
+        if (Type(value) == "AssociativeList") {
             models.codon.MG_REV.set_branch_length.p = value;
             if (models.codon.MG_REV.set_branch_length.p / terms.branch_length) {
                 models.codon.MG_REV.set_branch_length.p[terms.branch_length] =  models.codon.MG_REV.set_branch_length.p[terms.branch_length] * 3;
@@ -131,20 +131,20 @@ function models.codon.MG_REV.set_branch_length(model, value, parameter) {
         } else {
             models.codon.MG_REV.set_branch_length.p = 3*value;
         }
-    
+
         return models.generic.SetBranchLength(model, models.codon.MG_REV.set_branch_length.p, parameter);
     }
 
 
     models.codon.MG_REV.set_branch_length.beta = model.generic.GetLocalParameter(model, terms.parameters.nonsynonymous_rate);
     models.codon.MG_REV.set_branch_length.alpha = model.generic.GetLocalParameter(model, terms.parameters.synonymous_rate);
-    
+
     if ( null != models.codon.MG_REV.set_branch_length.beta &&  null != models.codon.MG_REV.set_branch_length.alpha) {
-    
+
 
         models.codon.MG_REV.set_branch_length.alpha.p = parameter + "." + models.codon.MG_REV.set_branch_length.alpha;
         models.codon.MG_REV.set_branch_length.beta.p = parameter + "." + models.codon.MG_REV.set_branch_length.beta;
-    
+
         if (Type(value) == "AssociativeList") {
             if (value[terms.model.branch_length_scaler] == terms.model.branch_length_constrain) {
                 if (parameters.IsIndependent(models.codon.MG_REV.set_branch_length.alpha.p)) {
@@ -171,14 +171,14 @@ 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;
 
                 }
-                                
+
                 if (parameters.IsIndependent(models.codon.MG_REV.set_branch_length.beta.p)) {
                     Eval(models.codon.MG_REV.set_branch_length.beta.p + ":=(" + value[terms.model.branch_length_scaler] + ")*" + value[terms.branch_length]);
                     models.codon.MG_REV.set_branch_length.lp += 1;
@@ -187,13 +187,13 @@ function models.codon.MG_REV.set_branch_length(model, value, parameter) {
             }
         } 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.alpha.p)) { // alpha 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, "");
-                
-                
+
+
                     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);
                     parameters.RemoveConstraint(models.codon.MG_REV.set_branch_length.beta);
@@ -201,21 +201,21 @@ function models.codon.MG_REV.set_branch_length(model, value, parameter) {
                 } else { // beta IS constrained
                     //parameters.SetConstraint (models.codon.MG_REV.set_branch_length.beta, parameters.GetConstraint (models.codon.MG_REV.set_branch_length.alpha.p),"");
                     //parameters.SetConstraint (models.codon.MG_REV.set_branch_length.alpha, models.codon.MG_REV.set_branch_length.alpha.p,"");
-                
-                    /** the branch length expression is going to be in terms of ** template ** parameters, but constraints will be in 
+
+                    /** the branch length expression is going to be in terms of ** template ** parameters, but constraints will be in
                         terms of instantiated parameters, so the expression to solve for needs to be temporarily bound to the global variable **/
 
-                
                     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, "");
-                
+
                     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);
                     parameters.RemoveConstraint ( models.codon.MG_REV.set_branch_length.beta);
                     parameters.RemoveConstraint ( models.codon.MG_REV.set_branch_length.alpha.p);
-                                    
+
+
                     messages.log ("models.codon.MG_REV.set_branch_length: " + models.codon.MG_REV.set_branch_length.alpha.p + "=" + models.codon.MG_REV.set_branch_length.lp);
-                
+                    model["models.codon.MG_REV.set_branch_length"] = {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/MG_REV_PROPERTIES.bf b/res/TemplateBatchFiles/libv3/models/codon/MG_REV_PROPERTIES.bf
index 2be76eacd..bf69ee216 100644
--- a/res/TemplateBatchFiles/libv3/models/codon/MG_REV_PROPERTIES.bf
+++ b/res/TemplateBatchFiles/libv3/models/codon/MG_REV_PROPERTIES.bf
@@ -1,5 +1,5 @@
 RequireVersion ("2.5.7");
- 
+
 LoadFunctionLibrary("../codon.bf");
 LoadFunctionLibrary("../DNA.bf");
 LoadFunctionLibrary("../parameters.bf");
@@ -19,10 +19,11 @@ LoadFunctionLibrary("../protein.bf");
                         'one letter amino-acid' -> property
                     }
  */
- 
 
- 
-terms.model.residue_name_map = "residue_to_id_map"; 
+
+
+terms.model.residue_name_map = "residue_to_id_map";
+terms.model.MG_REV_PROP.mean_prop = "mean_property_value";
 
 
 terms.model.MG_REV_PROP.predefined = {
@@ -140,123 +141,123 @@ terms.model.MG_REV_PROP.predefined = {
      },
      "LCAP" : {
         "Chemical Composition" : {
-            {-0.5868927454579467} 
-            {-0.5868460291700075} 
-            {-0.5868457923690357} 
-            {-0.586845591345447} 
-            {-0.5868454280017823} 
-            {1.442413127765356} 
-            {-0.02935869926260699} 
-            {0.4272890876194427} 
-            {-0.5868451756885411} 
-            {-0.3012025870591999} 
-            {0.2419262642781204} 
-            {0.6848994365857728} 
-            {1.312947120977898} 
-            {-0.1156419909884284} 
-            {1.384845017074097} 
-            {0.7270716263181436} 
-            {3.342731546150394} 
-            {-0.4010042423293907} 
-            {0.3415170202317384} 
-            {0.4699025925311998} 
+            {-0.5868927454579467}
+            {-0.5868460291700075}
+            {-0.5868457923690357}
+            {-0.586845591345447}
+            {-0.5868454280017823}
+            {1.442413127765356}
+            {-0.02935869926260699}
+            {0.4272890876194427}
+            {-0.5868451756885411}
+            {-0.3012025870591999}
+            {0.2419262642781204}
+            {0.6848994365857728}
+            {1.312947120977898}
+            {-0.1156419909884284}
+            {1.384845017074097}
+            {0.7270716263181436}
+            {3.342731546150394}
+            {-0.4010042423293907}
+            {0.3415170202317384}
+            {0.4699025925311998}
             },
         "Polarity" : {
-            {-0.9956290326072353} 
-            {-1.107220007612595} 
-            {-0.9956292005349376} 
-            {-0.8100842857525447} 
-            {-0.7357356892564315} 
-            {0.4912036016809371} 
-            {0.04436400722996305} 
-            {0.2672975016276415} 
-            {0.08234475255665361} 
-            {-0.624379765568966} 
-            {0.9376066737521076} 
-            {0.9747140546900138} 
-            {1.383229010106464} 
-            {1.271792102330754} 
-            {1.904064473242525} 
-            {1.643233506555707} 
-            {-0.8846717344020271} 
-            {-0.9219033329616851} 
-            {0.9747140551029528} 
-            {0.4163765565466709} 
+            {-0.9956290326072353}
+            {-1.107220007612595}
+            {-0.9956292005349376}
+            {-0.8100842857525447}
+            {-0.7357356892564315}
+            {0.4912036016809371}
+            {0.04436400722996305}
+            {0.2672975016276415}
+            {0.08234475255665361}
+            {-0.624379765568966}
+            {0.9376066737521076}
+            {0.9747140546900138}
+            {1.383229010106464}
+            {1.271792102330754}
+            {1.904064473242525}
+            {1.643233506555707}
+            {-0.8846717344020271}
+            {-0.9219033329616851}
+            {0.9747140551029528}
+            {0.4163765565466709}
             },
-            
+
         "Volume" : {
-            {1.257637282334762} 
-            {0.7687291690973896} 
-            {0.7687283894814801} 
-            {0.6290409911059951} 
-            {0.1401332510120709} 
-            {-1.070491916990394} 
-            {-1.058850814502372} 
-            {-0.3953354747308322} 
-            {-1.093771866791469} 
-            {1.350763381619076} 
-            {0.4195123000264476} 
-            {0.1634169621566906} 
-            {-0.5117403202399798} 
-            {0.9549794083529999} 
-            {-0.5583026469068318} 
-            {0.1168546775994032} 
-            {-0.535022037570798} 
-            {2.142318464323941} 
-            {1.071384899476407} 
-            {-1.745647766340383} 
+            {1.257637282334762}
+            {0.7687291690973896}
+            {0.7687283894814801}
+            {0.6290409911059951}
+            {0.1401332510120709}
+            {-1.070491916990394}
+            {-1.058850814502372}
+            {-0.3953354747308322}
+            {-1.093771866791469}
+            {1.350763381619076}
+            {0.4195123000264476}
+            {0.1634169621566906}
+            {-0.5117403202399798}
+            {0.9549794083529999}
+            {-0.5583026469068318}
+            {0.1168546775994032}
+            {-0.535022037570798}
+            {2.142318464323941}
+            {1.071384899476407}
+            {-1.745647766340383}
             },
         "Iso-electric Point" : {
-            {-0.3057542421356372} 
-            {-0.02273844767848476} 
-            {-9.737035422433548e-05} 
-            {-0.1585867248775605} 
-            {-0.03405832632980638} 
-            {-0.192548376206011} 
-            {0.1583919144982974} 
-            {0.07914812431496633} 
-            {-0.01141750903277056} 
-            {-0.2038687480243035} 
-            {0.8885757410301605} 
-            {-0.2095283729063603} 
-            {-0.3453753360155064} 
-            {2.105545575730344} 
-            {-1.839706056109631} 
-            {-1.584990785596673} 
-            {-0.5378282523391487} 
-            {-0.07368049231084244} 
-            {2.682902758933077} 
-            {-0.02839872447370546} 
+            {-0.3057542421356372}
+            {-0.02273844767848476}
+            {-9.737035422433548e-05}
+            {-0.1585867248775605}
+            {-0.03405832632980638}
+            {-0.192548376206011}
+            {0.1583919144982974}
+            {0.07914812431496633}
+            {-0.01141750903277056}
+            {-0.2038687480243035}
+            {0.8885757410301605}
+            {-0.2095283729063603}
+            {-0.3453753360155064}
+            {2.105545575730344}
+            {-1.839706056109631}
+            {-1.584990785596673}
+            {-0.5378282523391487}
+            {-0.07368049231084244}
+            {2.682902758933077}
+            {-0.02839872447370546}
             },
         "Hydropathy" : {
-            {1.301968439295016} 
-            {1.636753002222122} 
-            {1.871099553721135} 
-            {1.000656143286975} 
-            {1.770671640107803} 
-            {0.09673165440139944} 
-            {-0.171098473395166} 
-            {0.1302104025093898} 
-            {0.9671789981061145} 
-            {-0.0706624429784429} 
-            {-0.7067573940958518} 
-            {-0.8071939952622403} 
-            {-0.807193639057381} 
-            {-0.9411095669881511} 
-            {-0.8071943355554064} 
-            {-0.8071934015198535} 
-            {1.201528247390539} 
-            {0.06325356334516005} 
-            {-1.14198047362342} 
-            {0.2306464200526206} 
+            {1.301968439295016}
+            {1.636753002222122}
+            {1.871099553721135}
+            {1.000656143286975}
+            {1.770671640107803}
+            {0.09673165440139944}
+            {-0.171098473395166}
+            {0.1302104025093898}
+            {0.9671789981061145}
+            {-0.0706624429784429}
+            {-0.7067573940958518}
+            {-0.8071939952622403}
+            {-0.807193639057381}
+            {-0.9411095669881511}
+            {-0.8071943355554064}
+            {-0.8071934015198535}
+            {1.201528247390539}
+            {0.06325356334516005}
+            {-1.14198047362342}
+            {0.2306464200526206}
             }
      }
 };
 
- 
+
 lfunction models.codon.MG_REV_PROP.ModelDescription(type, code, properties) {
 
-     
+
     // piggyback on the standard MG_REV model for most of the code
 
     mg_base = models.codon.MG_REV.ModelDescription (type, code);
@@ -265,6 +266,7 @@ lfunction models.codon.MG_REV_PROP.ModelDescription(type, code, properties) {
     mg_base[utility.getGlobalValue("terms.model.residue_properties")] = models.codon.MG_REV_PROP._munge_properties(properties);
     mg_base[utility.getGlobalValue("terms.model.residue_name_map")] = parameters.ValidateIDs (utility.Keys (mg_base[utility.getGlobalValue("terms.model.residue_properties")]));
     mg_base[utility.getGlobalValue("terms.model.post_definition")] = "models.codon.MG_REV_PROP.post_definition";
+    mg_base[utility.getGlobalValue("terms.model.set_branch_length")] = "models.codon.MG_REV_PROP.set_branch_length";
 
     return mg_base;
 }
@@ -272,33 +274,33 @@ lfunction models.codon.MG_REV_PROP.ModelDescription(type, code, properties) {
 /**
  * @name models.codon.MG_REV_PROP._munge_properties
  * @param {String or AssociativeList} properties
- * @return {AssociativeList} 
+ * @return {AssociativeList}
  * Convert a set of amino-acid properties into the expected format
  */
 
 
 lfunction models.codon.MG_REV_PROP._munge_properties (properties) {
- 
+
     if (Type (properties) == "String") {
         assert (^"terms.model.MG_REV_PROP.predefined" / properties, properties + " is not a valid predefined property set");
         return models.codon.MG_REV_PROP._munge_properties ((^"terms.model.MG_REV_PROP.predefined")[properties]);
     }
 
     assert (Type (properties) == "AssociativeList", "Properties definition must be an AssociativeArray or a String for predefined properties");
-    
+
         /*
              each property must map to either a dictionary with '1 letter AA' -> value
              or to an array with 20 values
         */
-         
-         
+
+
     valid_aa = utility.MatrixToDict (^"models.protein.alphabet");
-    
+
     for (key, value; in; properties) {
         if (Type (value) == "AssociativeList") {
             io.CheckAssertion ("utility.Array1D(`&value`)==20", "A dictionary of amino-acid properties must have 20 entries for " + key);
             for (_i = 0; _i < 20; _i += 1) {
-                assert (value / (^"models.protein.alphabet")[_i], 
+                assert (value / (^"models.protein.alphabet")[_i],
                     "Property " + key + " is not defined for amino-acid " + (^"models.protein.alphabet")[_i]);
             }
         } else {
@@ -314,7 +316,7 @@ lfunction models.codon.MG_REV_PROP._munge_properties (properties) {
             }
         }
     }
-    
+
     assert (utility.Array1D(properties) > 0, "At least one valid property set must be defined");
     return properties;
 }
@@ -368,21 +370,21 @@ lfunction models.codon.MG_REV_PROP._GenerateRate_generic (fromChar, toChar, name
 
 
         rate_entry = nuc_rate;
-        
+
         if (_tt[fromChar] != _tt[toChar]) {
             if (model_type == utility.getGlobalValue("terms.global")) {
                 aa_rate = {};
                 prop_count = Abs (properties);
                 prop_names = utility.Keys (properties);
                 (_GenerateRate.p[model_type]) [^"terms.parameters.log_omega_ratio"] = parameters.ApplyNameSpace (omega, namespace);
-        
+
                 for (prop_name, prop_values; in; properties) {
                     prop_diff = Abs (prop_values[_tt[fromChar]] - prop_values[_tt[toChar]]);
-                    term_rate = parameters.ApplyNameSpace(lambda + "_" + property_id_map[prop_name], namespace); 
+                    term_rate = parameters.ApplyNameSpace(lambda + "_" + property_id_map[prop_name], namespace);
                     (_GenerateRate.p[model_type])[terms.propertyImportance (prop_name)] = term_rate;
                     aa_rate + ("`term_rate`*" + prop_diff);
                 }
-                
+
                 //aa_rate = parameters.ApplyNameSpace(lambda, namespace);
                 rate_entry += "*Min(10000,Exp(" + (_GenerateRate.p[model_type]) [^"terms.parameters.log_omega_ratio"] + "-(" + Join("+",aa_rate) + ")))";
              } else {
@@ -390,10 +392,10 @@ lfunction models.codon.MG_REV_PROP._GenerateRate_generic (fromChar, toChar, name
                 prop_count = Abs (properties);
                 prop_names = utility.Keys (properties);
                 (_GenerateRate.p[model_type]) [beta_term] = beta;
-        
+
                 for (prop_name, prop_values; in; properties) {
                     prop_diff = Abs (prop_values[_tt[fromChar]] - prop_values[_tt[toChar]]);
-                    term_rate = lambda + "_" + property_id_map[prop_name]; 
+                    term_rate = lambda + "_" + property_id_map[prop_name];
                     (_GenerateRate.p[model_type])[terms.propertyImportance (prop_name)] = term_rate;
                     aa_rate + ("`term_rate`*" + prop_diff);
                 }
@@ -425,4 +427,35 @@ lfunction  models.codon.MG_REV_PROP.post_definition (model) {
         parameters.SetRange(id, prop_range);
         parameters.SetValue(id, 0.1);
     }
+
+    models.generic.post.definition (model);
+}
+
+lfunction models.codon.MG_REV_PROP.set_branch_length(model, value, parameter) {
+
+    if (utility.Has (model, ^"terms.model.MG_REV_PROP.mean_prop", "Number")) {
+      properties = model.GetLocalParameters_RegExp(model, terms.propertyImportance (""));
+      for (tag, id; in; properties) {
+        parameters.SetValue (id, model[^"terms.model.MG_REV_PROP.mean_prop"]);
+      }
+      if (utility.Has (model, "fraction_same", "Number")) {
+        fs = model["fraction_same"];
+        models.codon.MG_REV.set_branch_length(model,value,parameter);
+        weighted = model ["models.codon.MG_REV.set_branch_length"];
+        if (Type (weighted) == "AssociativeList") {
+          for (tag, id; in; properties) {
+            parameters.SetValue (id, 0.);
+          }
+          models.codon.MG_REV.set_branch_length(model,value,parameter);
+          unweighted = model ["models.codon.MG_REV.set_branch_length"];
+          for (tag, value; in; weighted) {
+            //console.log (fs * unweighted[tag] + (1-fs) * value);
+            parameters.SetValue (tag, fs * unweighted[tag] + (1-fs) * value);
+          }
+          return 0;
+        }
+      }
+    }
+    return models.codon.MG_REV.set_branch_length(model,value,parameter);
+
 }
diff --git a/res/TemplateBatchFiles/libv3/models/model_functions.bf b/res/TemplateBatchFiles/libv3/models/model_functions.bf
index d1085a84d..c974d4305 100644
--- a/res/TemplateBatchFiles/libv3/models/model_functions.bf
+++ b/res/TemplateBatchFiles/libv3/models/model_functions.bf
@@ -26,6 +26,22 @@ lfunction model.GetParameters_RegExp(model, re) {
     return result;
 }
 
+lfunction model.GetLocalParameters_RegExp(model, re) {
+
+    names = utility.Filter (utility.Keys ((model[utility.getGlobalValue ("terms.parameters")])[ utility.getGlobalValue("terms.local")]),
+                            "_parameter_description_",
+                            "None != regexp.Find (_parameter_description_, `&re`)");
+
+    result = {};
+    count  = utility.Array1D (names);
+    for (k = 0; k < count; k += 1) {
+        result [names[k]] = ((model[utility.getGlobalValue ("terms.parameters")])[ utility.getGlobalValue("terms.local")])[names[k]];
+    }
+
+    return result;
+}
+
+
 /**
  * @name model.ApplyModelToTree
  * @param id {String}
@@ -43,7 +59,7 @@ function model.ApplyModelToTree (id, tree, model_list, rules) {
 	    // OR
 	    // DEFAULT : model id
 
-    
+
 	    if (Abs (rules["DEFAULT"])) {
             ExecuteCommands ("UseModel (" + rules["DEFAULT"] + ");
                               Tree `id` = " + tree["string"] + ";
@@ -90,13 +106,13 @@ function model.ApplyModelToTree (id, tree, model_list, rules) {
 
 	} else {
 	    // TO DO: REMOVE HARDCODING
-	    
-	    
+
+
 		model.ApplyModelToTree.modelID = model_list[model_list ["INDEXORDER"][0]];
 		ExecuteCommands ("UseModel (" + model.ApplyModelToTree.modelID[terms.id] + ");
 						  Tree `id` = " + tree["string"] + ";
 						  ");
-						  
+
 	}
 }
 
@@ -223,7 +239,7 @@ function model.generic.DefineModel (model_spec, id, arguments, data_filter, esti
 	model.generic.DefineModel.model = utility.CallFunction (model_spec, arguments);
 
 
- 
+
 	// Add data filter information to model description
 	if ( None != data_filter) {
 	    models.generic.AttachFilter (model.generic.DefineModel.model, data_filter);
@@ -234,8 +250,8 @@ function model.generic.DefineModel (model_spec, id, arguments, data_filter, esti
 
 
     // Define type of frequency estimator
-    
-    
+
+
 	if (None != estimator_type) {
 		model.generic.DefineModel.model [terms.model.frequency_estimator] = estimator_type;
 	}
@@ -351,7 +367,7 @@ function models.generic.ConstrainBranchLength (model, value, parameter) {
     if (Type (value) == "Number") {
         if (Abs((model[terms.parameters])[terms.local]) == 1) {
             if (Type (model [terms.model.branch_length_string]) == "String") {
-            
+
                 models.generic.ConstrainBranchLength.expression = model [terms.model.branch_length_string];
                 models.generic.ConstrainBranchLength.bl = (Columns ((model[terms.parameters])[terms.local]))[0];
                 models.generic.ConstrainBranchLength.bl.p = parameter + "." + models.generic.ConstrainBranchLength.bl;
@@ -382,8 +398,8 @@ function models.generic.ConstrainBranchLength (model, value, parameter) {
  * @returns the number of constraints generated (0 or 1)
  */
 function models.generic.SetBranchLength (model, value, parameter) {
-    
-    
+
+
      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];
@@ -537,7 +553,7 @@ lfunction models.BindGlobalParameters (models, filter) {
 
 
     if (Type (models) == "AssociativeList" && utility.Array1D (models) > 1) {
-    
+
         reference_set = (((models[0])[utility.getGlobalValue("terms.parameters")])[utility.getGlobalValue("terms.global")]);
         candidate_set = utility.UniqueValues(utility.Filter (utility.Keys (reference_set), "_key_",
             "regexp.Find (_key_,`&filter`)"
diff --git a/res/TemplateBatchFiles/libv3/models/parameters.bf b/res/TemplateBatchFiles/libv3/models/parameters.bf
index dc39d75d9..5de46b415 100644
--- a/res/TemplateBatchFiles/libv3/models/parameters.bf
+++ b/res/TemplateBatchFiles/libv3/models/parameters.bf
@@ -524,7 +524,6 @@ lfunction parameters.GetConstraint(parameter) {
 function parameters.SetConstraint(id, value, global_tag) {
     if (Type(id) == "String") {
         if (Abs(id)) {
-            //console.log ("`global_tag` `id` := " + value);
             ExecuteCommands("`global_tag` `id` := " + value);
         }
     } else {
diff --git a/res/TemplateBatchFiles/libv3/tasks/alignments.bf b/res/TemplateBatchFiles/libv3/tasks/alignments.bf
index 98f195e4f..ecf679625 100644
--- a/res/TemplateBatchFiles/libv3/tasks/alignments.bf
+++ b/res/TemplateBatchFiles/libv3/tasks/alignments.bf
@@ -103,7 +103,7 @@ lfunction alignments.ReadNucleotideDataSet_aux(dataset_name) {
             partitions = Transpose(dfpm);
         }
     }
-    
+
     return {
         utility.getGlobalValue("terms.data.sequences"): Eval("`dataset_name`.species"),
         utility.getGlobalValue("terms.data.sites"): Eval("`dataset_name`.sites"),
@@ -405,6 +405,22 @@ function alignments.ReadNucleotideAlignment(file_name, dataset_name, datafilter_
     return data_info;
 }
 
+/**
+ * Take an input filter and check if it has duplicate sequences
+ * @name alignments.HasDuplicateSequences
+ * @param {String} filter_in - The name of an existing filter
+ * @param {Number} check_mode -
+      -1 : exact match
+      -2 : exact match + gaps match non-gaps");
+      -3 : superset filtering
+      -4 : partial match filtering
+ - @returns the number of *duplicate* sequences
+ */
+lfunction alignments.HasDuplicateSequences (filter_in, check_mode) {
+  GetDataInfo (duplicate_info, ^filter_in, check_mode);
+  return duplicate_info["UNIQUE_SEQUENCES"];
+}
+
 /**
  * Take an input filter, replace all identical sequences with a single copy
  * Optionally, rename the sequences to indicate copy # by adding ':copies'
@@ -416,7 +432,7 @@ function alignments.ReadNucleotideAlignment(file_name, dataset_name, datafilter_
  */
 lfunction alignments.CompressDuplicateSequences (filter_in, filter_out, rename) {
 
-    GetDataInfo (duplicate_info, ^filter_in, -2);    
+    GetDataInfo (duplicate_info, ^filter_in, -2);
     DataSetFilter ^filter_out = CreateFilter (^filter_in, 1, "", Join (",", duplicate_info["UNIQUE_INDICES"]));
 
     if (rename) {
@@ -516,6 +532,7 @@ lfunction alignments.TranslateCodonsToAminoAcids (sequence, offset, code) {
     l = Abs (sequence);
 	translation = "";
 	translation * (l/3+1);
+
 	for (k = offset; k < l; k += 3) {
 		codon = sequence[k][k+2];
 		if (code [utility.getGlobalValue("terms.code.mapping")] / codon) {
@@ -525,7 +542,7 @@ lfunction alignments.TranslateCodonsToAminoAcids (sequence, offset, code) {
 		    if (codon == "---") {
 			    translation * "-";
 		    } else {
-			    translation * "?";
+		        translation * "?";
 			}
 	    }
 	}
@@ -953,9 +970,9 @@ lfunction alignment.MapCodonsToAA(codon_sequence, aa_sequence, this_many_mm, map
 
 lfunction alignment.ExportPartitionedNEXUS (filter, breakPoints, trees, file, isCodon) {
     utility.ToggleEnvVariable ("DATA_FILE_PRINT_FORMAT", 4);
-    
+
     fprintf (file, CLEAR_FILE, KEEP_OPEN, ^filter, "\n");
-    
+
     breakPointsCount = utility.Array1D (breakPoints);
     partCount      = breakPointsCount + 1;
     currentStart    = 0;
@@ -974,10 +991,10 @@ lfunction alignment.ExportPartitionedNEXUS (filter, breakPoints, trees, file, is
             }
             currentEnd = currentEnd - 1;
         }
- 
+
         fprintf (file, "\tCHARSET span_", p + 1, " = ", currentStart + 1, "-", currentEnd + 1, ";\n");
-        
- 
+
+
         if (!lastPartition) {
             currentStart = breakPoints[p] + 1;
         }
@@ -987,7 +1004,7 @@ lfunction alignment.ExportPartitionedNEXUS (filter, breakPoints, trees, file, is
     for (p = 0; p < partCount; p += 1) {
         fprintf (file, "\tTREE tree_", p + 1, " = ", trees[p], ";\n");
     }
-    
+
     fprintf (file, "END;\n");
     utility.ToggleEnvVariable ("DATA_FILE_PRINT_FORMAT", None);
 }
diff --git a/res/TemplateBatchFiles/libv3/tasks/mpi.bf b/res/TemplateBatchFiles/libv3/tasks/mpi.bf
index 9360f345e..b15da82d8 100644
--- a/res/TemplateBatchFiles/libv3/tasks/mpi.bf
+++ b/res/TemplateBatchFiles/libv3/tasks/mpi.bf
@@ -314,6 +314,7 @@ namespace mpi {
                                            "1" : jobs [i],
                                            "2" : &scores}, callback);
         }
+        
 
         Call (callback, -1, Call (handler, lf_id, jobs[0], &scores), {"0" : lf_id, "1" : jobs [0], "2" : &scores});
 
@@ -339,7 +340,6 @@ namespace mpi {
         task_ids = utility.Keys (tasks);
         task_count = Abs (tasks);
         for (i = 0; i < task_count; i+=1) {
-
             parameters.SetValues (tasks[task_ids[i]]);
             LFCompute (^lf_id, ll);
             results [task_ids[i]] = ll;
diff --git a/res/TemplateBatchFiles/libv3/tasks/trees.bf b/res/TemplateBatchFiles/libv3/tasks/trees.bf
index 9ee63d0e4..59fdc1b6c 100644
--- a/res/TemplateBatchFiles/libv3/tasks/trees.bf
+++ b/res/TemplateBatchFiles/libv3/tasks/trees.bf
@@ -102,7 +102,7 @@ lfunction trees.GetTreeString._sanitize(string) {
     if (utility.GetEnvVariable("_KEEP_I_LABELS_")) {
         utility.ToggleEnvVariable("INTERNAL_NODE_PREFIX", None);
     }
-    
+
     return string;
 }
 
@@ -135,9 +135,9 @@ lfunction trees.GetTreeString(look_for_newick_tree) {
             utility.SetEnvVariable ("LAST_FILE_IO_EXCEPTION", None);
             utility.ToggleEnvVariable ("SOFT_FILE_IO_EXCEPTIONS",TRUE);
             fscanf(PROMPT_FOR_FILE, REWIND, "Raw", treeString);
-            
+
             look_for_newick_tree = utility.getGlobalValue ("LAST_FILE_PATH");
-            
+
             if (None != utility.GetEnvVariable ("LAST_FILE_IO_EXCEPTION")) {
                 if (utility.getGlobalValue ("LAST_RAW_FILE_PROMPT") ==  utility.getGlobalValue ("terms.trees.neighbor_joining")) {
                     datafilter_name = utility.GetEnvVariable(utility.getGlobalValue ("terms.trees.data_for_neighbor_joining"));
@@ -148,11 +148,11 @@ lfunction trees.GetTreeString(look_for_newick_tree) {
                 }
                 utility.SetEnvVariable ("LAST_FILE_IO_EXCEPTION", None);
             }
-            
+
             utility.ToggleEnvVariable ("SOFT_FILE_IO_EXCEPTIONS",None);
             fprintf(stdout, "\n");
 
-  
+
             if (regexp.Find(treeString, "^#NEXUS")) {
                 ExecuteCommands(treeString);
 
@@ -254,7 +254,13 @@ lfunction trees.LoadAnnotatedTopology(look_for_newick_tree) {
 lfunction trees.LoadAnnotatedTopologyAndMap(look_for_newick_tree, mapping) {
 
     reverse = {};
-    utility.ForEach(utility.Keys(mapping), "_key_", "`&reverse`[`&mapping`[_key_]] = _key_");
+
+    for (k,v; in; mapping) {
+      reverse[v] = k;
+    }
+
+    //utility.ForEach(utility.Keys(mapping), "_key_", "`&reverse`[`&mapping`[_key_]] = _key_");
+    
     io.CheckAssertion("Abs (`&mapping`) == Abs (`&reverse`)", "The mapping between original and normalized tree sequence names must be one to one");
     utility.ToggleEnvVariable("TREE_NODE_NAME_MAPPING", reverse);
     result = trees.ExtractTreeInfo(trees.GetTreeString(look_for_newick_tree));
@@ -354,7 +360,7 @@ lfunction trees.branch_names(tree, respect_case) {
 /**
  * @name trees.KillZeroBranches
  * Given a tree dict and a list of matching parameter estimates
- * this returns a modified tree dict with all zero-branch length 
+ * this returns a modified tree dict with all zero-branch length
  * internal branches removed and modifies in in place
  * @param tree - dict of the tree object
  * @param estimates - dict with branch length estimates
@@ -362,9 +368,9 @@ lfunction trees.branch_names(tree, respect_case) {
  * @param zero_internal -- store branches that were deleted here
  * @return modified tree
  */
- 
+
 lfunction trees.KillZeroBranches (tree, estimates, branch_set, zero_internal) {
-    
+
     for (branch, value; in; tree [^"terms.trees.partitioned"]) {
         if (value == ^"terms.tree_attributes.internal") {
             if (estimates / branch) {
@@ -384,11 +390,11 @@ lfunction trees.KillZeroBranches (tree, estimates, branch_set, zero_internal) {
             }
         }
         return trees.ExtractTreeInfoFromTopology (&T);
-        
+
     }
-    
+
     return tree;
-    
+
 }
 
 /**
@@ -431,8 +437,8 @@ lfunction trees.RootTree(tree_info, root_on) {
  * - list of models
  */
 lfunction trees.ExtractTreeInfoFromTopology(topology_object) {
-    
- 
+
+
     branch_lengths = BranchLength(^topology_object, -1);
     branch_names   = BranchName(^topology_object, -1);
     branch_count   = utility.Array1D (branch_names) - 1;
@@ -500,7 +506,7 @@ lfunction trees.ExtractTreeInfo(tree_string) {
 
     Topology T = tree_string;
     return trees.ExtractTreeInfoFromTopology (&T);
-    
+
 }
 
 /**
@@ -993,11 +999,11 @@ lfunction tree._GenerateRandomTreeDraw2 (nodes) {
     do {
         n2 = Random (0,n) $ 1;
     } while (n1 == n2);
-    
-    
+
+
     nodes - r[n1];
     nodes - r[n2];
-    
+
     return {"0" : +r[n1], "1" : +r[n2]};
 }
 
@@ -1011,18 +1017,18 @@ lfunction tree.GenerateRandomTree (N, rooted, branch_name, branch_length) {
 
     total_nodes = N + internal_nodes;
     flat_tree  = {total_nodes, 4}["-1"];
-    
+
 
     available_to_join = {};
     for (k = 0; k < N; k+=1) {
         available_to_join[k] = TRUE;
     }
-    
+
 
     current_parent_node = N;
     downto = 1 + (rooted == 0);
-    
-    
+
+
     while (Abs (available_to_join) > downto) {
         pair = tree._GenerateRandomTreeDraw2 (available_to_join);
         flat_tree[pair[0]][0] = current_parent_node;
@@ -1032,13 +1038,13 @@ lfunction tree.GenerateRandomTree (N, rooted, branch_name, branch_length) {
         available_to_join[current_parent_node] = TRUE;
         current_parent_node += 1;
     }
-    
+
     if (!rooted) { // attach the last node to the root
         available_to_join - (total_nodes-1);
         leaf_index = +((Rows(available_to_join))[0]);
         flat_tree[leaf_index][0] = total_nodes-1;
         flat_tree[total_nodes-1][3] = leaf_index;
-        
+
     }
 
     return tree._NewickFromMatrix (&flat_tree, total_nodes-1, branch_name, branch_length);
@@ -1068,7 +1074,7 @@ lfunction tree._NewickFromMatrix (flat_tree, index, branch_name, branch_length)
             } else {
                 bn := "";
             }
-            
+
             return "("   +  tree._NewickFromMatrix (flat_tree, (^flat_tree)[index][1], branch_name, branch_length) +
                      "," +  tree._NewickFromMatrix (flat_tree, (^flat_tree)[index][2], branch_name, branch_length) +
                      ")" + bn + bl;
diff --git a/src/core/batchlan.cpp b/src/core/batchlan.cpp
index 1136e332c..b80eca26f 100644
--- a/src/core/batchlan.cpp
+++ b/src/core/batchlan.cpp
@@ -3191,6 +3191,12 @@ void      _ElementaryCommand::ExecuteCase52 (_ExecutionList& chain) {
         SetStatusLine ("Simulating Data");
         { // lf must be deleted before the referenced datafilters
             _LikelihoodFunction lf (filter_specification, nil);
+            
+            /*_SimpleList gl;
+            lf.GetGlobalVars(gl);
+            gl.Each ([] (long vi, unsigned long) -> void { StringToConsole(*LocateVar(vi)->GetName()); NLToConsole();});
+            */
+            
             lf.Simulate (*sim_dataset, exclusions, category_values, category_names, root_states, do_internals?(main_file?&spool_file:&kEmptyString):nil);
             SetStatusLine ("Idle");
         }
@@ -4137,7 +4143,8 @@ bool       _ElementaryCommand::MakeGeneralizedLoop  (_String*p1, _String*p2, _St
             // None != iterator_value, but we don't know
             
             _StringBuffer iterator_condition;
-            iterator_condition << "None!=" << (_String*)advance->parameters.GetItem (advance->parameters.countitems() -1);
+            iterator_condition << kEndIteration.Enquote() << "!=" << (_String*)advance->parameters.GetItem (advance->parameters.countitems() -1);
+            //StringToConsole(iterator_condition);NLToConsole();
             target.GetIthCommand(for_return)->MakeJumpCommand (&iterator_condition, for_return+1, target.lLength,target);
             
         } else {
diff --git a/src/core/batchlanruntime.cpp b/src/core/batchlanruntime.cpp
index 875fe13f8..edcfc8ee2 100644
--- a/src/core/batchlanruntime.cpp
+++ b/src/core/batchlanruntime.cpp
@@ -1873,7 +1873,7 @@ bool      _ElementaryCommand::HandleAdvanceIterator(_ExecutionList& current_prog
            }
           
            if (row >= source_object->GetHDim()) {
-               ((_Variable*)parameters.GetItem (reciever_count << 1))->SetValue (new _MathObject, false, false, NULL);
+               ((_Variable*)parameters.GetItem (reciever_count << 1))->SetValue (new _FString (kEndIteration), false, false, NULL); // end loop
            } else {
                ((_Variable*)parameters.GetItem (reciever_count << 1))->SetValue (source_object->GetMatrixCell(row, column), false, false, NULL);
                if (reciever_count == 2) {
@@ -1907,7 +1907,7 @@ bool      _ElementaryCommand::HandleAdvanceIterator(_ExecutionList& current_prog
                       ((_Variable*)parameters.GetItem (reciever_count+1))->SetValue ((HBLObjectRef)state.get_object(), false, false,NULL);
                   }
               } else {
-                  ((_Variable*)parameters.GetItem (reciever_count << 1))->SetValue (new _MathObject, false, false,NULL);
+                  ((_Variable*)parameters.GetItem (reciever_count << 1))->SetValue (new _FString (kEndIteration), false, false,NULL);
                   // iterator done
               }
           } else {
@@ -1935,7 +1935,7 @@ bool      _ElementaryCommand::HandleAdvanceIterator(_ExecutionList& current_prog
                   }
                   simpleParameters[2]++;
               } else {
-                ((_Variable*)parameters.GetItem (reciever_count << 1))->SetValue (new _MathObject, false, false,NULL);
+                ((_Variable*)parameters.GetItem (reciever_count << 1))->SetValue (new _FString (kEndIteration), false, false,NULL);
               }
               
           }
@@ -3471,16 +3471,26 @@ bool      _ElementaryCommand::HandleGetString (_ExecutionList& current_program)
               }
               else {  // formula string
 
-                if (index1 == -3) {
+                if (index1 == -3 || index1 == -4) {
                   _StringBuffer local, global;
+                    
                   _SimpleList var_index;
-                  var_index << var->get_index ();
-                  if (var->IsIndependent()) {
-                      //printf ("ExportIndVariables\n");
-                    ExportIndVariables (global, local, &var_index);
+                  if (index1 == -3 || var->IsIndependent())  {
+                      var_index << var->get_index ();
+                      if (var->IsIndependent()) {
+                          //printf ("ExportIndVariables\n");
+                        ExportIndVariables (global, local, &var_index);
+                      } else {
+                          //printf ("ExportDepVariables\n");
+                        ExportDepVariables(global, local, &var_index);
+                      }
                   } else {
-                      //printf ("ExportDepVariables\n");
-                    ExportDepVariables(global, local, &var_index);
+                      _AVLList vl (&var_index);
+                      var->ScanForVariables(vl, true);
+                      _SimpleList ind_vars = var_index.Filter([] (long index, unsigned long) -> bool {return LocateVar(index)->IsIndependent();}),
+                                  dep_vars = var_index.Filter([] (long index, unsigned long) -> bool {return !LocateVar(index)->IsIndependent();});
+                      ExportIndVariables(global, local, &ind_vars);
+                      ExportDepVariables(global, local, &dep_vars);
                   }
                   return_value = make_fstring_pointer ( & ((*new _StringBuffer (128L)) << global << local << '\n'));
 
diff --git a/src/core/calcnode.cpp b/src/core/calcnode.cpp
index 1f40387b8..39d3455d6 100644
--- a/src/core/calcnode.cpp
+++ b/src/core/calcnode.cpp
@@ -54,6 +54,9 @@
 using namespace hy_global;
 using namespace hy_env;
 
+//#define _UBER_VERBOSE_MX_UPDATE_DUMP
+//#define _UBER_VERBOSE_MX_UPDATE_DUMP_EVAL 1
+
 //_______________________________________________________________________________________________
 
 _CalcNode::_CalcNode    () {
@@ -480,12 +483,12 @@ bool        _CalcNode::RecomputeMatrix  (long categID, long totalCategs, _Matrix
     }
   
     #ifdef _UBER_VERBOSE_MX_UPDATE_DUMP
-      if (1|| likeFuncEvalCallCount == _UBER_VERBOSE_MX_UPDATE_DUMP_LF_EVAL && gVariables) {
+     // if (1|| likeFuncEvalCallCount == _UBER_VERBOSE_MX_UPDATE_DUMP_LF_EVAL && gVariables) {
         for (unsigned long i=0; i<gVariables->lLength; i++) {
           _Variable* curVar = LocateVar(gVariables->GetElement(i));
-          fprintf (stderr, "[_CalcNode::RecomputeMatrix] Node %s, var %s, value = %15.12g\n", GetName()->sData, curVar->GetName()->sData, curVar->Compute()->Value());
+          fprintf (stderr, "[_CalcNode::RecomputeMatrix] Node %s, var %s, value = %15.12g\n", GetName()->get_str(), curVar->GetName()->get_str(), curVar->Compute()->Value());
         }
-      }
+      //}
     #endif
 
     /*
@@ -501,7 +504,7 @@ bool        _CalcNode::RecomputeMatrix  (long categID, long totalCategs, _Matrix
     if (!storeRateMatrix) {
       if (totalCategs>1) {
   #ifdef _UBER_VERBOSE_MX_UPDATE_DUMP
-        fprintf (stderr, "[_CalcNode::RecomputeMatrix] Deleting category %ld for node %s at %p\n", categID, GetName()->sData, GetCompExp(categID));
+        fprintf (stderr, "[_CalcNode::RecomputeMatrix] Deleting category %ld for node %s at %p\n", categID, GetName()->get_str(), GetCompExp(categID));
   #endif
         if (clear_exponentials()) {
             DeleteObject(GetCompExp(categID, true));
@@ -520,7 +523,7 @@ bool        _CalcNode::RecomputeMatrix  (long categID, long totalCategs, _Matrix
     if (isExplicitForm && bufferedOps) {
         _Matrix * bufferedExp = (_Matrix*)GetExplicitFormModel()->Compute (0,nil, bufferedOps);
         #ifdef _UBER_VERBOSE_MX_UPDATE_DUMP
-            fprintf (stderr, "[_CalcNode::RecomputeMatrix] Setting (buffered) category %ld/%ld for node %s\n", categID, totalCategs, GetName()->sData);
+            fprintf (stderr, "[_CalcNode::RecomputeMatrix] Setting (buffered) category %ld/%ld for node %s\n", categID, totalCategs, GetName()->get_str());
          #endif
         SetCompExp ((_Matrix*)bufferedExp->makeDynamic(), totalCategs>1?categID:-1);
         return false;
@@ -575,7 +578,7 @@ bool        _CalcNode::RecomputeMatrix  (long categID, long totalCategs, _Matrix
             }
 
             #ifdef _UBER_VERBOSE_MX_UPDATE_DUMP
-                fprintf (stderr, "[_CalcNode::RecomputeMatrix] Setting category %ld/%ld for node %s\n", categID, totalCategs, GetName()->sData);
+                fprintf (stderr, "[_CalcNode::RecomputeMatrix] Setting category %ld/%ld for node %s\n", categID, totalCategs, GetName()->get_str());
             #endif
             SetCompExp ((_Matrix*)(isExplicitForm?temp:temp->Exponentiate(1., true)), totalCategs>1?categID:-1);
 
diff --git a/src/core/dataset.cpp b/src/core/dataset.cpp
index 159f3df2d..6f0aa12dd 100644
--- a/src/core/dataset.cpp
+++ b/src/core/dataset.cpp
@@ -515,6 +515,7 @@ void _DataSet::Finalize(void) {
 
       for (long i3 = 0; i3 < lLength; i3++) {
         tC = (_Site *)(*(_List *)this)(i3);
+        tC->TrimSpace();
         tC->SetRefNo(0);
       }
       if (dsh) {
@@ -1862,19 +1863,20 @@ void ReadNextLine (FILE* fp, _String *s, FileState* fs, bool, bool upCase) {
     char lastc;
     
     if (fp) {
-        lastc = fgetc(fp);
+        lastc = getc_unlocked (fp);
     } else {
         lastc = fs->pInSrc<fs->theSource->length()?fs->theSource->char_at(fs->pInSrc++):0;
     }
     
+    
     if (fs->fileType != 3) { // not NEXUS - do not skip [..]
         if (fp)
-            while ( !feof(fp) && lastc!=10 && lastc!=13 ) {
+            while ( !feof_unlocked(fp) && lastc!=10 && lastc!=13 ) {
                 if (lastc) {
                     tempBuffer << lastc;
                 }
                 
-                lastc = fgetc(fp);
+                lastc = getc_unlocked(fp);
             }
         else
             while (lastc && lastc!=10 && lastc!=13 ) {
@@ -1887,7 +1889,7 @@ void ReadNextLine (FILE* fp, _String *s, FileState* fs, bool, bool upCase) {
             lastc = toupper(lastc);
         }
         
-        while (((fp&&(!feof(fp)))||(fs->theSource&&(fs->pInSrc<=fs->theSource->length ()))) && lastc!='\r' && lastc!='\n') {
+        while (((fp&&(!feof_unlocked(fp)))||(fs->theSource&&(fs->pInSrc<=fs->theSource->length ()))) && lastc!='\r' && lastc!='\n') {
             if (lastc=='[') {
                 if (fs->isSkippingInNEXUS) {
                     ReportWarning ("Nested comments in NEXUS really shouldn't be used.");
@@ -1908,7 +1910,7 @@ void ReadNextLine (FILE* fp, _String *s, FileState* fs, bool, bool upCase) {
                 if (upCase) {
                     lastc = toupper(fgetc(fp));
                 } else {
-                    lastc = fgetc(fp);
+                    lastc = getc_unlocked(fp);
                 }
             } else {
                 if (upCase) {
@@ -1927,7 +1929,7 @@ void ReadNextLine (FILE* fp, _String *s, FileState* fs, bool, bool upCase) {
     
     tempBuffer.TrimSpace();
     
-    if ( (fp && feof(fp)) || (fs->theSource && fs->pInSrc >= fs->theSource->length()) ) {
+    if ( (fp && feof_unlocked (fp)) || (fs->theSource && fs->pInSrc >= fs->theSource->length()) ) {
         if (tempBuffer.empty ()) {
             *s = "";
             return;
@@ -1972,11 +1974,18 @@ _DataSet* ReadDataSetFile (FILE*f, char execBF, _String* theS, _String* bfName,
     
     try {
     
-        _String         CurrentLine = hy_env::data_file_tree_string & "={{}};",
+        if (f) flockfile (f);
+        
+        hy_env::EnvVariableSet(hy_env::data_file_tree_string, new _Matrix, false);
+        
+        _String savedLine;
+        
+        /*_String         CurrentLine = hy_env::data_file_tree_string & "={{}};",
                         savedLine;
         
         _ExecutionList reset (CurrentLine);
-        reset.Execute();
+        reset.Execute();*/
+        
     #ifdef __HYPHYMPI__
         if (hy_mpi_node_rank == 0L)
     #endif
@@ -2030,12 +2039,12 @@ _DataSet* ReadDataSetFile (FILE*f, char execBF, _String* theS, _String* bfName,
         }
     #endif
         
-        
+        _String     CurrentLine;
         
         //if (f==NULL) return (_DataSet*)result.makeDynamic();
         // nothing to do
         
-        CurrentLine = kEmptyString;
+        //CurrentLine = kEmptyString;
         
         ReadNextLine (f,&CurrentLine,&fState);
         if (CurrentLine.empty()) {
@@ -2278,6 +2287,7 @@ _DataSet* ReadDataSetFile (FILE*f, char execBF, _String* theS, _String* bfName,
             {
                 _TranslationTable *trialTable = new _TranslationTable (hy_default_translation_table);
                 trialTable->baseLength = 2;
+                if (f) funlockfile (f);
                 _DataSet * res2 = ReadDataSetFile (f, execBF, theS, bfName, namespaceID, trialTable);
                 if (res2->GetNoTypes()) {
                     DeleteObject (result);
@@ -2334,10 +2344,11 @@ _DataSet* ReadDataSetFile (FILE*f, char execBF, _String* theS, _String* bfName,
         }
     } catch (const _String& err) {
         DeleteObject (result);
+        if (f) funlockfile (f);
         HandleApplicationError(err);
         result = nil;
     }
-    
+    if (f) funlockfile (f);
     return result;
 }
 
diff --git a/src/core/dataset_filter.cpp b/src/core/dataset_filter.cpp
index 07f4add5d..a479ffe23 100644
--- a/src/core/dataset_filter.cpp
+++ b/src/core/dataset_filter.cpp
@@ -108,73 +108,69 @@ unsigned long    _DataSetFilter::FindUniqueSequences  (_SimpleList& indices, _Si
     seqs   = theNodeMap.lLength,
     unit   = GetUnitLength();
     
+    _SimpleList _sequenceHashes;
+    _AVLListXL  sequenceHashes     (&_sequenceHashes);
+    
+    auto      check_sequence_match = [&] (long f, long rawSequenceIdx, _SimpleList* & sequencesWithSameScore) -> long {
+        if (f >= 0) {
+            sequencesWithSameScore = (_SimpleList*)sequenceHashes.GetXtra (f);
+            for (long k = 0; k<sequencesWithSameScore->lLength; k++) {
+                bool fit = true;
+                f = sequencesWithSameScore->list_data[k];
+                
+                long fRaw = theNodeMap.list_data[indices.list_data[f]];
+                
+                for (unsigned long site = 0; site < sites && fit; site++){
+                    _Site * thisSite = theData->GetSite(site);
+                    if (thisSite->char_at(fRaw) != thisSite->char_at(rawSequenceIdx)) {
+                        fit = false;
+                        break;
+                    }
+                }
+                
+                if (fit) {
+                    map << f;
+                    counts.list_data[f] ++;
+                    return f;
+                }
+            }
+        }
+        return -1;
+    };
+    
+    auto      add_sequence_match = [&] (_SimpleList* sameScore, long sequenceHash, long sequenceIndex, bool do_map)->void {
+        if (!sameScore) {
+            sameScore = new _SimpleList;
+            sequenceHashes.Insert ((BaseRef)sequenceHash,(long)sameScore,false);
+        }
+        
+        (*sameScore) << indices.lLength;
+        
+        if (do_map) {
+            map     << indices.lLength;
+            indices << sequenceIndex;
+            counts  << 1;
+        }
+    };
+
     if (match_mode == kUniqueMatchExact) { // exact match
-        _SimpleList hashSupport;
-        _AVLListXL  sequenceHashes     (&hashSupport);
         
         for (unsigned long sequenceIndex = 0; sequenceIndex < seqs; sequenceIndex ++){
             _String * thisSequence = GetSequenceCharacters (sequenceIndex);
-            
+  
             long     sequenceHash   = thisSequence->Adler32(),
-            f              = sequenceHashes.Find ((BaseRef)sequenceHash),
+            f              = sequenceHashes.FindLong (sequenceHash),
             rawSequenceIdx = theNodeMap.list_data[sequenceIndex];
             
             DeleteObject (thisSequence);
-            
+
             _SimpleList * sameScore = nil;
-            if (f>=0) {
-                sameScore = (_SimpleList*)sequenceHashes.GetXtra (f);
-                for (long k = 0; k<sameScore->lLength; k++) {
-                    bool fit = true;
-                    f = sameScore->list_data[k];
-                    
-                    long fRaw = theNodeMap.list_data[indices.list_data[f]];
-                    
-                    for (unsigned long site = 0; site < sites && fit; site++){
-                        for (unsigned long unitIndex = 0; unitIndex < unit; unitIndex ++){
-                            _Site * thisSite = theData->GetSite(theMap.list_data[unit*site+unitIndex]);
-                            if (thisSite->char_at(fRaw) != thisSite->char_at(rawSequenceIdx)) {
-                                fit = false;
-                                break;
-                            }
-                        }
-                    }
-                    
-                    if (fit) {
-                        /*
-                        StringToConsole (*GetSequenceName(sequenceIndex));
-                        BufferToConsole("==");
-                        StringToConsole (*GetSequenceName(f));
-                        NLToConsole();
-                        
-                        StringToConsole (*GetSequenceCharacters(sequenceIndex));
-                        NLToConsole();
-                        StringToConsole (*GetSequenceCharacters(f));
-                        NLToConsole();
-                        NLToConsole();
-                        */
-
-                        map << f;
-                        counts.list_data[f] ++;
-                        
-                    } else {
-                        f = -1;
-                    }
-                }
-            }
+            f = check_sequence_match (f, rawSequenceIdx, sameScore);
+            
             if (f==-1) { // fit failed or unique site
-                if (!sameScore) {
-                    sameScore = new _SimpleList;
-                    sequenceHashes.Insert ((BaseRef)sequenceHash,(long)sameScore,false);
-                }
-                
-                (*sameScore) << indices.lLength;
-                map     << indices.lLength;
-                indices << sequenceIndex;
-                counts  << 1;
+                add_sequence_match (sameScore, sequenceHash, sequenceIndex, true);
             }
         }
-        
     }
     else{
         long             vd  = GetDimension(true);
@@ -182,13 +178,37 @@ unsigned long    _DataSetFilter::FindUniqueSequences  (_SimpleList& indices, _Si
         hyFloat      *translatedVector  = new hyFloat [vd],
                      *translatedVector2 = new hyFloat [vd];
         
-        _String         state1 (unit), state2 (unit);
+        _String      state1 (unit), state2 (unit);
+        
+        long         *referenceStates = (long*) new long[sites];
         
         sites = sites / unit;
         
         for (long sequenceIndex = 0; sequenceIndex < seqs; sequenceIndex++) {
             bool check_state = false;
             
+            //printf ("%ld (%ld)\n", sequenceIndex, indices.countitems());
+            
+            _String * thisSequence = GetSequenceCharacters (sequenceIndex);
+            
+            long     sequenceHash   = thisSequence->Adler32(),
+                     f              = sequenceHashes.FindLong (sequenceHash),
+                     rawSequenceIdx = theNodeMap.list_data[sequenceIndex];
+            
+            DeleteObject(thisSequence);
+            _SimpleList * sameScore = nil;
+            f = check_sequence_match (f, rawSequenceIdx, sameScore);
+            
+            if (f >= 0L) {
+                //printf ("Adler Hash match for %ld (%ld/%ld)\n", sequenceIndex, f, indices.countitems());
+                continue;
+            }
+            
+            for (long m=0L; m<sites; m++) {
+                RetrieveState (m,sequenceIndex,    state1,false);
+                referenceStates[m] = Translate2Frequencies (state1, translatedVector,  true);
+            }
+            
             for (long idx=0L; idx<indices.countitems(); idx++) {
                 check_state = true;
                 for (long m=0L; m<sites; m++) {
@@ -196,100 +216,92 @@ unsigned long    _DataSetFilter::FindUniqueSequences  (_SimpleList& indices, _Si
                     RetrieveState (m,sequenceIndex,    state1,false);
                     RetrieveState (m,indices.get(idx), state2,false);
                     
-                    long idx1 = Translate2Frequencies (state1, translatedVector,  true),
-                         idx2 = Translate2Frequencies (state2, translatedVector2, true);
-                    
-                    /*if (idx == 85 && sequenceIndex == 88) {
-                        printf ("(%ld, %ld) %ld = %ld %ld\n", sequenceIndex, indices.list_data[idx], m, idx1, idx2);
-                    }*/
-                    
-                    if (idx2 >=0 && idx1 >=0) { // fully resolved
-                        if (idx1==idx2) {
-                            continue;
-                        } else {
-                            check_state = false;
-                            break;
-                        }
-                    } else {
-                        
-                        // check for equal ambigs
-                        long k = 0L;
-                        for (; k < vd; k++){
-                            if (translatedVector[k] != translatedVector2[k]){
-                                break;
-                            }
-                        }
-                        
-                        if (k == vd) { // resolutions matched -- this is always OK
-                            continue;
-                        }
+                    if (state1 != state2) { // only check different states
+                        long idx1 = referenceStates[m],//Translate2Frequencies (state1, translatedVector,  true),
+                             idx2 = Translate2Frequencies (state2, translatedVector2, true);
                         
+                        /*if (idx == 85 && sequenceIndex == 88) {
+                            printf ("(%ld, %ld) %ld = %ld %ld\n", sequenceIndex, indices.list_data[idx], m, idx1, idx2);
+                        }*/
                         
-                        if (match_mode == kUniqueMatchExactOrGap) {
-                            
-                            long count1 = 0L,
-                                 count2 = 0L;
-                            
-                            for (long t = 0L; t<vd; t++) {
-                                count1 += translatedVector[t]  > 0.0;
-                                count2 += translatedVector2[t] > 0.0;
-                            }
-                            
-                            // Gaps can match resolved positions, but not the other way around
-                            // count2 is the number of resolutions in the currently stored sequence
-                            if (count2 > 1L && count2 < vd || count1 < vd) {
+                        if (idx2 >=0 && idx1 >=0) { // fully resolved
+                            if (idx1==idx2) {
+                                continue;
+                            } else {
                                 check_state = false;
                                 break;
                             }
+                        } else {
+                            //RetrieveState (m,sequenceIndex,    state1,false);
+                            Translate2Frequencies (state1, translatedVector,  true);
+                            // check for equal ambigs
+                            long k = 0L;
+                            for (; k < vd; k++){
+                                if (translatedVector[k] != translatedVector2[k]){
+                                    break;
+                                }
+                            }
                             
+                            if (k == vd) { // resolutions matched -- this is always OK
+                                continue;
+                            }
                             
-                        } else {
-                            bool first  = (match_mode == kUniqueMatchSuperset),
-                                 second = first;
                             
-                            if (match_mode == kUniqueMatchSuperset){
-                                for (long t = 0L; (first||second)&&(t<vd); t++) {
-                                    if (translatedVector[t]>0.0) {
-                                        second &= (translatedVector2[t]>0.0);
-                                    }
-                                    if (translatedVector2[t]>0.0) {
-                                        first  &= (translatedVector[t]>0.0);
-                                    }
+                            if (match_mode == kUniqueMatchExactOrGap) {
+                                
+                                long count1 = 0L,
+                                     count2 = 0L;
+                                
+                                for (long t = 0L; t<vd; t++) {
+                                    count1 += translatedVector[t]  > 0.0;
+                                    count2 += translatedVector2[t] > 0.0;
                                 }
-                                if (!(first||second)) {
+                                
+                                // Gaps can match resolved positions, but not the other way around
+                                // count2 is the number of resolutions in the currently stored sequence
+                                if (count2 > 1L && count2 < vd || count1 < vd) {
                                     check_state = false;
                                     break;
                                 }
+                                
+                                
                             } else {
-                                for (long t = 0L; t<vd; t++) {
-                                    if (translatedVector[t]>0.0) {
-                                        second |= (translatedVector2[t]>0.0);
+                                bool first  = (match_mode == kUniqueMatchSuperset),
+                                     second = first;
+                                
+                                if (match_mode == kUniqueMatchSuperset){
+                                    for (long t = 0L; (first||second)&&(t<vd); t++) {
+                                        if (translatedVector[t]>0.0) {
+                                            second &= (translatedVector2[t]>0.0);
+                                        }
+                                        if (translatedVector2[t]>0.0) {
+                                            first  &= (translatedVector[t]>0.0);
+                                        }
                                     }
-                                    if (translatedVector2[t]>0.0) {
-                                        first  |= (translatedVector[t]>0.0);
+                                    if (!(first||second)) {
+                                        check_state = false;
+                                        break;
+                                    }
+                                } else {
+                                    for (long t = 0L; t<vd; t++) {
+                                        if (translatedVector[t]>0.0) {
+                                            second |= (translatedVector2[t]>0.0);
+                                        }
+                                        if (translatedVector2[t]>0.0) {
+                                            first  |= (translatedVector[t]>0.0);
+                                        }
+                                    }
+                                    if (!(first&&second)) {
+                                        check_state = false;
+                                        break;
                                     }
-                                }
-                                if (!(first&&second)) {
-                                    check_state = false;
-                                    break;
                                 }
                             }
                         }
                     }
                 }
-                
+            
                 if (check_state) {
-                    /*StringToConsole (*GetSequenceName(sequenceIndex));
-                    BufferToConsole("== [ref]");
-                    StringToConsole (*GetSequenceName(indices.get (idx)));
-                    NLToConsole();
-                    
-                    StringToConsole (*GetSequenceCharacters(sequenceIndex));
-                    NLToConsole();
-                    StringToConsole (*GetSequenceCharacters(indices.get (idx)));
-                    NLToConsole();
-                    NLToConsole();*/
-
                     map << idx;
                     counts.list_data[idx] ++;
                     break;
@@ -297,6 +309,7 @@ unsigned long    _DataSetFilter::FindUniqueSequences  (_SimpleList& indices, _Si
             }
             
             if (!check_state){
+                add_sequence_match (sameScore, sequenceHash, sequenceIndex, false);
                 map     << indices.lLength;
                 indices << sequenceIndex;
                 counts  << 1;
@@ -306,6 +319,7 @@ unsigned long    _DataSetFilter::FindUniqueSequences  (_SimpleList& indices, _Si
         
         delete [] translatedVector;
         delete [] translatedVector2;
+        delete [] referenceStates;
     }
     
     
@@ -812,17 +826,6 @@ _String* _DataSetFilter::MakeSiteBuffer (void) const {
     return new _String ((unsigned long)unitLength);
 }
 
-//_______________________________________________________________________
-void _DataSetFilter::retrieve_individual_site_from_raw_coordinates (_String& store, unsigned long site, unsigned long sequence) const {
-  if (unitLength==1UL) {
-    store.set_char (0, (((_String**)theData->list_data)[theData->theMap.list_data[theMap.list_data[site]]])->char_at (sequence));
-  } else {
-    site*=unitLength;
-    for (unsigned long k = 0UL; k<unitLength; k++) {
-      store.set_char (k, ((_String**)theData->list_data)[theData->theMap.list_data[theMap.list_data[site++]]]->char_at(sequence));
-    }
-  }
-}
 
 //_______________________________________________________________________
 
@@ -868,11 +871,6 @@ _List *  _DataSetFilter::ComputePatternToSiteMap (void) const {
 }
 
 
-//_______________________________________________________________________
-
-char _DataSetFilter::direct_index_character (unsigned long site, unsigned long sequence) const {
-  return (((_String**)theData->list_data)[theData->theMap.list_data[theMap.list_data[site]]])->char_at(sequence);
-}
 
 //_______________________________________________________________________
  
@@ -1021,14 +1019,19 @@ long    _DataSetFilter::HasExclusions (unsigned long site, _SimpleList* theExc,
     if (theNodeMap.countitems()) {
         _String buffer ((unsigned long)GetUnitLength());
   
+        unsigned long   filter_dim = GetDimension(false);
+
         for (unsigned long k = 0UL; k<theNodeMap.countitems(); k++) {
             RetrieveState           (site, k, buffer, false);
-            Translate2Frequencies   (buffer, store, false, false);
+            long idx =              Translate2Frequencies   (buffer, store, false, false);
             
-            unsigned long   filter_dim = GetDimension(false);
             long            found_forbidden = -1;
-          
-            for (unsigned long character = 0UL ;  character < filter_dim;  character ++) {
+            if (idx >= 0L) {
+                if (theExc->Find(idx) != kNotFound) {
+                    found_forbidden = k;
+                }
+            } else {
+                for (unsigned long character = 0UL ;  character < filter_dim;  character ++) {
                    if (store[character] > 0.0) {
                       if (theExc->Find(character) == kNotFound) {
                           found_forbidden = -1;
@@ -1037,6 +1040,7 @@ long    _DataSetFilter::HasExclusions (unsigned long site, _SimpleList* theExc,
                         found_forbidden = k;
                       }
                   }
+                }
             }
           
             if (found_forbidden >= 0L) {
@@ -1621,13 +1625,16 @@ long    _DataSetFilter::Translate2Frequencies (_String const& str, hyFloat* parv
     long mapped_resolution_count = correct_for_exclusions && theExclusions.nonempty() ? theExclusions.CorrectForExclusions(store, resolution_count) : resolution_count;
     
     /* handle the cases when no unambiguous resolutions were available */
+
+    if (mapped_resolution_count == 1L) {
+        parvect[store[0]] = 1.;
+        return store[0];
+    }
+    
     for (long i = 0L; i < mapped_resolution_count; i++) {
         parvect[store[i]] = 1.;
     }
     
-    if (mapped_resolution_count == 1L) {
-        return store[0];
-    }
     if (mapped_resolution_count == 0L && resolution_count == 0L && smear) {
         InitializeArray(parvect, dimension, 1.);
     }
diff --git a/src/core/formula.cpp b/src/core/formula.cpp
index 4e8670434..8deae2f31 100644
--- a/src/core/formula.cpp
+++ b/src/core/formula.cpp
@@ -1842,8 +1842,14 @@ long      _Formula::ExtractMatrixExpArguments (_List* storage) {
                       * next_op = GetIthTerm(i + 1UL);
 
               if (! cacheUpdated && next_op->CanResultsBeCached(this_op)) {
+                  /*
+                   StringToConsole("\n----\n");
+                  ObjectToConsole(FetchVar(LocateVarByName("k"))->Compute());
+                  NLToConsole();*/
+
                   _Stack temp;
                   this_op->Execute (temp);
+                  
 
                   _Matrix *currentArg  = (_Matrix*)temp.Pop(true),
                           *cachedArg   = (_Matrix*)((HBLObjectRef)(*resultCache)(cacheID)),
@@ -2827,7 +2833,7 @@ void    _Formula::ConvertToTree (bool err_msg) {
             } else {
                 theTree = (node<long>*)nodeStack(0);
             }
-        } catch (_String const e) {
+        } catch (_String const& e) {
             if (err_msg) {
                 HandleApplicationError(e);
             }
diff --git a/src/core/global_things.cpp b/src/core/global_things.cpp
index 158401ba6..c0afb193e 100644
--- a/src/core/global_things.cpp
+++ b/src/core/global_things.cpp
@@ -120,11 +120,12 @@ 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.15"),
+                     kHyPhyVersion  = _String ("2.5.16"),
     
                     kNoneToken = "None",
                     kNullToken = "null",
-                    kNoKWMatch = "__input_value_not_given__";
+                    kNoKWMatch = "__input_value_not_given__",
+                    kEndIteration = "__iterator_end__loop__";
   
     _String
                      hy_base_directory,
diff --git a/src/core/include/dataset_filter.h b/src/core/include/dataset_filter.h
index 8de150571..0412296ec 100644
--- a/src/core/include/dataset_filter.h
+++ b/src/core/include/dataset_filter.h
@@ -352,8 +352,19 @@ class _DataSetFilter : public BaseObj {
   void internalToStr(FILE *, _StringBuffer *);
   
   
-   void retrieve_individual_site_from_raw_coordinates (_String & store_here, unsigned long site, unsigned long sequence) const;
-   inline char direct_index_character (unsigned long site, unsigned long sequence) const;
+   inline void retrieve_individual_site_from_raw_coordinates (_String & store, unsigned long site, unsigned long sequence) const {
+        if (unitLength==1UL) {
+          store.set_char (0, (((_String**)theData->list_data)[theData->theMap.list_data[theMap.list_data[site]]])->char_at (sequence));
+        } else {
+          site*=unitLength;
+          for (unsigned long k = 0UL; k<unitLength; k++) {
+            store.set_char (k, ((_String**)theData->list_data)[theData->theMap.list_data[theMap.list_data[site++]]]->char_at(sequence));
+          }
+        }
+   }
+    inline char direct_index_character (unsigned long site, unsigned long sequence) const {
+        return (((_String**)theData->list_data)[theData->theMap.list_data[theMap.list_data[site]]])->char_at(sequence);
+    }
   
    _String *accessCache;
 
diff --git a/src/core/include/function_templates.h b/src/core/include/function_templates.h
index b525a4b80..b83cdf5dc 100644
--- a/src/core/include/function_templates.h
+++ b/src/core/include/function_templates.h
@@ -156,6 +156,8 @@ void ArrayForEach(ARG_TYPE *array, unsigned long dimension,
 template <typename ARG_TYPE>
 void InitializeArray(ARG_TYPE *array, unsigned long dimension,
                      ARG_TYPE &&value) {
+  #pragma clang loop unroll_count(8)
+  #pragma GCC unroll 4
   for (unsigned long i = 0UL; i < dimension; i++) {
     array[i] = value;
   }
diff --git a/src/core/include/global_things.h b/src/core/include/global_things.h
index e505c05ac..826a82180 100644
--- a/src/core/include/global_things.h
+++ b/src/core/include/global_things.h
@@ -308,7 +308,8 @@ namespace hy_global {
   kErrorStringUnterminatedMatrix,
   kNoneToken,
   kNullToken,
-  kNoKWMatch;
+  kNoKWMatch,
+  kEndIteration;
   
   extern  _String  hy_base_directory,
   hy_error_log_name,
diff --git a/src/core/include/polynoml.h b/src/core/include/polynoml.h
index 3caf8b104..5afab462a 100644
--- a/src/core/include/polynoml.h
+++ b/src/core/include/polynoml.h
@@ -124,7 +124,8 @@ class _Polynomial : public _MathObject {
 
     _Polynomial             (void);
     _Polynomial             (_SimpleList&);
-    _Polynomial             (_Polynomial&);
+    _Polynomial             (_Polynomial const&);
+    _Polynomial const &     operator = (_Polynomial const&);
     _Polynomial             (hyFloat);
     _Polynomial             (_Variable&);
     virtual                 ~_Polynomial ();
diff --git a/src/core/include/tree.h b/src/core/include/tree.h
index ea85839d5..6cd88d4e8 100644
--- a/src/core/include/tree.h
+++ b/src/core/include/tree.h
@@ -132,7 +132,7 @@ class _TheTree: public _TreeTopology
     _List*      SampleAncestors                 (_DataSetFilter*, node<long>*);
     void        PurgeTree                       (void);
 
-    long        ComputeReleafingCost            (_DataSetFilter const*, long, long, _SimpleList* = nil, long = 0) const;
+    long        ComputeReleafingCost            (_DataSetFilter const*, long, long, _SimpleList* = nil, long = 0, _SimpleList const* = nil) const;
     long        ComputeReleafingCostChar        (_DataSetFilter const*, long, long, _SimpleList const*) const;
     void        DumpingOrder                    (_DataSetFilter*, _SimpleList&);
     void        SetTreeCodeBase                 (long);
diff --git a/src/core/likefunc.cpp b/src/core/likefunc.cpp
index 6cd4305f3..555404741 100644
--- a/src/core/likefunc.cpp
+++ b/src/core/likefunc.cpp
@@ -2309,7 +2309,7 @@ hyFloat  _LikelihoodFunction::Compute        (void)
         if (result >= 0.) {
             if (result >= __DBL_EPSILON__ * 1.e4) {
                 char buffer [2048];
-                snprintf (buffer, 2047, "Internal error: Encountered a positive log-likelihood (%g) at evaluation %ld, mode %ld, template %ld. This is usually a consequence of 'infinite-like' parameter values", result, likeFuncEvalCallCount-1, computeMode, templateKind);
+                snprintf (buffer, 2047, "Internal error: Encountered a positive log-likelihood (%g) at evaluation %ld, mode %d, template %ld. This is usually a consequence of 'infinite-like' parameter values", result, likeFuncEvalCallCount-1, computeMode, templateKind);
                 _TerminateAndDump(buffer, true);
             } else {
                 result = 0.;
@@ -5290,7 +5290,7 @@ long    _LikelihoodFunction::Bracket (long index, hyFloat& left, hyFloat& middle
 
     while (1) {
 
-        while (middle-leftStep < lowerBound) {
+        while (middle-leftStep <= lowerBound) {
             if (verbosity_level > 100) {
               char buf [512];
               snprintf (buf, sizeof(buf), "\n\t[_LikelihoodFunction::Bracket (index %ld) HANDLING LEFT BOUNDARY CASES] : LB = %g, current try = %.16g, current evaluated midpoint value = %.16g (%s)", index, lowerBound, middle-leftStep, middleValue, first ? "first" : "NOT first");
@@ -5337,7 +5337,7 @@ long    _LikelihoodFunction::Bracket (long index, hyFloat& left, hyFloat& middle
         }
 
 
-        while (rightStep+middle > upperBound) {
+        while (rightStep+middle >= upperBound) {
             rightStep = rightStep*.125; //MIN (rightStep*.125,upperBound - middle);
             if (rightStep<initialStep*.1 && index >= 0 || index < 0 && rightStep < STD_GRAD_STEP) {
                 if (!first) {
@@ -8135,6 +8135,9 @@ void    _LikelihoodFunction::Setup (bool check_reversibility)
                         alreadyDone = alreadyDoneModels.GetXtra (alreadyDone);
                     } else {
                         alreadyDone = IsModelReversible (treeModels.list_data[m]);
+                        if (!alreadyDone) {
+                            ReportWarning (_String ("Model ") & GetObjectNameByType(HY_BL_MODEL, treeModels.list_data[m])->Enquote() & " was determined to be non-reversible");
+                        }
                         alreadyDoneModels.Insert ((BaseRef)treeModels.list_data[m], alreadyDone);
                     }
                     if (alreadyDone) {
@@ -8148,6 +8151,10 @@ void    _LikelihoodFunction::Setup (bool check_reversibility)
                             _Matrix * my_freqs = (_Matrix*)f->GetValue();
                             if (!my_freqs->Equal (base_frequencies)) {
                                 isReversiblePartition = false;
+                                if (!alreadyDone) {
+                                    ReportWarning (_String ("Model ") & GetObjectNameByType(HY_BL_MODEL, treeModels.list_data[m])->Enquote() & " has a different base frequency vector compared to other models");
+                                }
+
                             }
                         }
                     }
@@ -8398,6 +8405,7 @@ hyFloat  _LikelihoodFunction::ComputeBlock (long index, hyFloat* siteRes, long c
               hyFloat correction = 0.;
                 
               for (blockID = 0; blockID < np; blockID ++)  {
+                //printf ("EVAL %ld BLOCK %ld (sites per %ld), RESULT %g\n", likeFuncEvalCallCount, blockID, sitesPerP, thread_results[blockID]);
                 if (thread_results[blockID] == -INFINITY) {
                   sum = -INFINITY;
                   break;
@@ -8720,6 +8728,10 @@ void    setComputingArrays (node<long>* startingNode, node<long>* childNode, _Si
 
 void        _LikelihoodFunction::OptimalOrder    (long index, _SimpleList& sl, _SimpleList const * clone) {
 
+    //printf ("\nEntered _LikelihoodFunction::OptimalOrder\n");
+    //TimeDifference tracker;
+    //tracker.Start();
+    
     _DataSetFilter const* df = GetIthFilter (index);
     _TheTree            * t   = GetIthTree(index);
     long                vLevel       = VerbosityLevel();
@@ -8769,7 +8781,7 @@ void        _LikelihoodFunction::OptimalOrder    (long index, _SimpleList& sl, _
                       child_count (t->get_flat_nodes().MapList([] (long n, unsigned long ) -> long {
                           return ((node <long>*)n)->get_num_nodes();
                       }));
-
+        
 
         while (completedSites<totalSites) {
             if (totalSites-completedSites<partition) {
@@ -8789,7 +8801,7 @@ void        _LikelihoodFunction::OptimalOrder    (long index, _SimpleList& sl, _
             } else {
                 for (k=completedSites+1; k<completedSites+partition; k++) {
                     partitionSites<<k;
-                    distances<<t->ComputeReleafingCost (df, completedSites, k);
+                    distances<<t->ComputeReleafingCost (df, completedSites, k, nil, 0, &child_count);
                     edges<<completedSites;
                 }
             }
@@ -8995,6 +9007,9 @@ void        _LikelihoodFunction::OptimalOrder    (long index, _SimpleList& sl, _
             }
         }
     }
+    
+    //printf ("Finished in %g seconds\n%s\n", tracker.TimeSinceStart (), buffer);
+
 
 }
 
@@ -9070,8 +9085,12 @@ void _LikelihoodFunction::SerializeLF(_StringBuffer & rec, char opt,
         }
     }
 
-    _String *lfName = (_String *)likeFuncNamesList(
+    _String const *lfName = (_String const *)likeFuncNamesList.GetItemRangeCheck(
         likeFuncList._SimpleList::Find((long) this));
+    
+    if (!lfName) { // inline LF
+        lfName = &kEmptyString;
+    }
 
     ReportWarning(_String("Serializing ") & *lfName & " with mode " &
                   _String((long) opt));
@@ -9115,6 +9134,10 @@ void _LikelihoodFunction::SerializeLF(_StringBuffer & rec, char opt,
             (long)(GetDataFilter(redirector->get(idx))
                        ->GetData()));
 
+        if (tIdx < 0) {
+            HandleApplicationError(_String("Can't serialize 'temporary' likelihood functions"));
+            return;
+        }
 
         tIdx = taggedDS.Insert((BaseRef) tIdx, taggedDS.countitems());
 
diff --git a/src/core/likefunc2.cpp b/src/core/likefunc2.cpp
index 24fa431df..748a91008 100644
--- a/src/core/likefunc2.cpp
+++ b/src/core/likefunc2.cpp
@@ -252,7 +252,7 @@ void            _LikelihoodFunction::SetupCategoryCaches      (void)
 
                 categoryTraversalTemplate.AppendNewInstance(container);
             }
-            catch (const _String error) {
+            catch (const _String& error) {
                 BatchDelete (catVarReferences,catVarCounts,catVarOffsets,hmmAndCOP,varType,container);
                 HandleApplicationError (error);
                 return;
diff --git a/src/core/matrix.cpp b/src/core/matrix.cpp
index 35206d976..91881d6ae 100644
--- a/src/core/matrix.cpp
+++ b/src/core/matrix.cpp
@@ -1907,130 +1907,137 @@ bool    _Matrix::IsValidTransitionMatrix() const {
 //_____________________________________________________________________________________________
 
 bool    _Matrix::IsReversible(_Matrix* freqs) {
-    if (!is_square() || (freqs && freqs->GetHDim () * freqs->GetVDim () != GetHDim())
-            || (!is_numeric() && !is_expression_based() ) ||
-            (freqs && !freqs->is_numeric() && !freqs->is_expression_based())) {
-        return false;
-    }
+    
+    try {
+        
+        if (!is_square()) {
+            throw _String ("Not a square matrix in _Matrix::IsReversible");
+        }
+        
+        if (freqs && freqs->GetHDim () * freqs->GetVDim () != GetHDim()) {
+            throw _String ("Incompatible frequency and rate matrix dimensions in _Matrix::IsReversible");
+        }
+        
+        if (!is_numeric() && !is_expression_based()) {
+            throw _String ("Unsupported rate matrix type in _Matrix::IsReversible");
+        }
 
-    bool   needAnalytics = is_expression_based() || (freqs && freqs->is_expression_based());
-    if (needAnalytics) {
-        if (freqs) {
-            for (long r = 0; r < hDim; r++)
-                for (long c = r+1; c < hDim; c++) {
-                    bool compResult = true;
-                    if (is_expression_based()) {
-                        _Formula* rc = GetFormula(r,c),
-                                  * cr = GetFormula(c,r);
+        if (freqs && !freqs->is_numeric() && !freqs->is_expression_based()) {
+            throw _String ("Unsupported frequency matrix type in _Matrix::IsReversible");
+        }
 
-                        if (rc && cr) {
-                            _Polynomial *rcp = (_Polynomial *)rc->ConstructPolynomial(),
-                                         *crp = (_Polynomial *)cr->ConstructPolynomial();
 
-                            if (rcp && crp) {
-                                HBLObjectRef     tr = nil,
-                                              tc = nil;
-
-                                if (freqs->is_expression_based()) {
-                                    if (freqs->GetFormula(r,0)) {
-                                        tr = freqs->GetFormula(r,0)->ConstructPolynomial();
-                                        if (tr) {
-                                            tr->AddAReference();
-                                        } else {
-                                            return false;
+        bool   needAnalytics = is_expression_based() || (freqs && freqs->is_expression_based());
+        if (needAnalytics) {
+            if (freqs) {
+                for (long r = 0; r < hDim; r++)
+                    for (long c = r+1; c < hDim; c++) {
+                        bool compResult = true;
+                        if (is_expression_based()) {
+                            _Formula* rc = GetFormula(r,c),
+                                      * cr = GetFormula(c,r);
+
+                            if (rc && cr) {
+                                _Polynomial *rcp = (_Polynomial *)rc->ConstructPolynomial(),
+                                             *crp = (_Polynomial *)cr->ConstructPolynomial();
+
+                                if (rcp && crp) {
+                                    HBLObjectRef     tr = nil,
+                                                     tc = nil;
+
+                                    if (freqs->is_expression_based()) {
+                                        if (freqs->GetFormula(r,0)) {
+                                            tr = freqs->GetFormula(r,0)->ConstructPolynomial();
+                                            if (tr) {
+                                                tr->AddAReference();
+                                            } else {
+                                                throw _String ("Could not convert matrix cell (") & r & ',' & c & ") to a polynomial";
+                                            }
                                         }
-                                    }
-                                    if (freqs->GetFormula(c,0)) {
-                                        tc = freqs->GetFormula(c,0)->ConstructPolynomial();
-                                        if (tc) {
-                                            tc->AddAReference();
-                                        } else {
-                                            DeleteObject (tr);
-                                            return false;
+                                        if (freqs->GetFormula(c,0)) {
+                                            tc = freqs->GetFormula(c,0)->ConstructPolynomial();
+                                            if (tc) {
+                                                tc->AddAReference();
+                                            } else {
+                                                DeleteObject (tr);
+                                                throw _String ("Could not convert frequency cell (") & c & ") to a polynomial";
+                                             }
                                         }
+                                    } else {
+                                        tr = new _Constant ((*freqs)[r]);
+                                        tc = new _Constant ((*freqs)[c]);
+                                    }
+                                    if (tr && tc) {
+                                        _Polynomial        * rcpF = (_Polynomial*)rcp->Mult(tr),
+                                                           * crpF = (_Polynomial*)crp->Mult(tc);
+
+                                        compResult         = rcpF->Equal(crpF);
+                                        DeleteObject (rcpF);
+                                        DeleteObject (crpF);
+                                        
+                                    } else {
+                                        compResult = !(tr||tc);
                                     }
-
-                                } else {
-                                    tr = new _Constant ((*freqs)[r]);
-                                    tc = new _Constant ((*freqs)[c]);
-                                }
-                                if (tr && tc) {
-                                    _Polynomial        * rcpF = (_Polynomial*)rcp->Mult(tr),
-                                                       * crpF = (_Polynomial*)crp->Mult(tc);
-
-                                    compResult         = rcpF->Equal(crpF);
                                     
-                                    /*if (!compResult) {
-                                        ObjectToConsole(rcpF);
-                                        NLToConsole();
-                                        ObjectToConsole(crpF);
-                                        NLToConsole();
-                                    }*/
+                                    DeleteObject (tr);
+                                    DeleteObject (tc);
                                     
-                                    DeleteObject (rcpF);
-                                    DeleteObject (crpF);
+                                    if (!compResult) {
+                                         throw _String ("Unequal cells at (") & r & ',' & c & ")";
+                                    }
                                 } else {
-                                    compResult = !(tr||tc);
+                                    throw _String ("Could not convert a matrix cell at (") & r & ',' & c & ") to a polynomial: " & _StringBuffer ((_StringBuffer*)rc->toStr(kFormulaStringConversionNormal));
                                 }
+                             } else {
+                                compResult = !(rc || cr);
+                            }
+                        }
+                    }
+            } else {
+                for (long r = 0; r < hDim; r++)
+                    for (long c = r+1; c < hDim; c++) {
+                        bool compResult = true;
+                        _Formula* rc = GetFormula(r,c),
+                                  * cr = GetFormula(c,r);
+
+                        if (rc && cr) {
+                            _Polynomial *rcp = (_Polynomial *)rc->ConstructPolynomial(),
+                                         *crp = (_Polynomial *)cr->ConstructPolynomial();
 
-                                DeleteObject (tr);
-                                DeleteObject (tc);
+                            if (rcp && crp) {
+                                compResult = rcp->Equal(crp);
                             } else {
-                                compResult = false;
+                                compResult = rc->EqualFormula(cr);
                             }
-
-                            //DeleteObject (rcp); DeleteObject (crp);
                         } else {
                             compResult = !(rc || cr);
                         }
-                    }
-                    if (!compResult) {
-                        return false;
-                    }
-                }
-        } else {
-            for (long r = 0; r < hDim; r++)
-                for (long c = r+1; c < hDim; c++) {
-                    bool compResult = true;
-                    _Formula* rc = GetFormula(r,c),
-                              * cr = GetFormula(c,r);
-
-                    if (rc && cr) {
-                        _Polynomial *rcp = (_Polynomial *)rc->ConstructPolynomial(),
-                                     *crp = (_Polynomial *)cr->ConstructPolynomial();
-
-                        if (rcp && crp) {
-                            compResult = rcp->Equal(crp);
-                        } else {
-                            compResult = rc->EqualFormula(cr);
-                        }
-
-                        //DeleteObject (rcp); DeleteObject (crp);
-                    } else {
-                        compResult = !(rc || cr);
-                    }
 
-                    if (!compResult) {
-                        return false;
-                    }
-                }
-        }
-        return true;
-    } else {
-        if (freqs) {
-            for (long r = 0; r < hDim; r++)
-                for (long c = r+1; c < hDim; c++)
-                    if (! CheckEqual ((*this)(r,c)*(*freqs)[r], (*this)(c,r)*(*freqs)[c])) {
-                        return false;
+                        if (!compResult) {
+                            return false;
+                        }
                     }
+            }
+            return true;
         } else {
-            for (long r = 0; r < hDim; r++)
-                for (long c = r+1; c < hDim; c++)
-                    if (! CheckEqual ((*this)(r,c), (*this)(c,r))) {
-                        return false;
-                    }
+            if (freqs) {
+                for (long r = 0; r < hDim; r++)
+                    for (long c = r+1; c < hDim; c++)
+                        if (! CheckEqual ((*this)(r,c)*(*freqs)[r], (*this)(c,r)*(*freqs)[c])) {
+                            return false;
+                        }
+            } else {
+                for (long r = 0; r < hDim; r++)
+                    for (long c = r+1; c < hDim; c++)
+                        if (! CheckEqual ((*this)(r,c), (*this)(c,r))) {
+                            return false;
+                        }
+            }
+            return true;
         }
-        return true;
+    } catch (const _String reason) {
+        ReportWarning (_String ("Reversibility checks failed: ") & reason);
+        return false;
     }
     return false;
 }
@@ -8337,7 +8344,7 @@ HBLObjectRef   _Matrix::MultinomialSample (_Constant *replicates) {
             }
         }
     }
-    catch (_String const err) {
+    catch (_String const& err) {
         HandleApplicationError (err);
     }
     return new _Matrix;
diff --git a/src/core/parser.cpp b/src/core/parser.cpp
index b003498f1..391dc9ecf 100644
--- a/src/core/parser.cpp
+++ b/src/core/parser.cpp
@@ -688,7 +688,7 @@ _Variable* CheckReceptacleCommandIDException (_String const* name, const long id
     }
   }
   
-  if (f<0L) {
+  {
     _Variable (*name, isGlobal);
     f = LocateVarByName (*name);
   }
diff --git a/src/core/polynoml.cpp b/src/core/polynoml.cpp
index d78550930..8ce1f1234 100644
--- a/src/core/polynoml.cpp
+++ b/src/core/polynoml.cpp
@@ -779,21 +779,28 @@ _Polynomial::_Polynomial (_Variable& v)
 
 //__________________________________________________________________________________
 
-_Polynomial::_Polynomial (_Polynomial& p)
-{
-    variableIndex.Duplicate (&p.variableIndex);
-    theTerms = new _PolynomialData;
-    if (p.theTerms) {
-        theTerms->Duplicate (p.theTerms);
-    } else {
-        theTerms->numberVars = variableIndex.countitems();
-    }
-    compList1.Duplicate (&p.compList1);
-    compList2.Duplicate (&p.compList2);
+_Polynomial::_Polynomial (_Polynomial const& p) {
+    *this = p;
 }
 
 //__________________________________________________________________________________
 
+_Polynomial const & _Polynomial::operator = (_Polynomial const& p) {
+    if (this != &p) {
+        variableIndex.Duplicate (&p.variableIndex);
+        theTerms = new _PolynomialData;
+        if (p.theTerms) {
+            theTerms->Duplicate (p.theTerms);
+        } else {
+            theTerms->numberVars = variableIndex.countitems();
+        }
+        compList1.Duplicate (&p.compList1);
+        compList2.Duplicate (&p.compList2);
+    }
+    return *this;
+}
+//__________________________________________________________________________________
+
 _Polynomial::~_Polynomial ()
 {
     if (theTerms) {
diff --git a/src/core/topology.cpp b/src/core/topology.cpp
index 7742a7b2a..3ed232b72 100644
--- a/src/core/topology.cpp
+++ b/src/core/topology.cpp
@@ -293,6 +293,8 @@ _AssociativeList*    _TreeTopology::MainTreeConstructor  (_String const& parms,
     
     static _String       kBootstrap ("bootstrap"),
                          kComment   ("comment");
+    
+    static long          kBLLookahead = 128;
   
     _SimpleList nodeStack,
     nodeNumbers;
@@ -477,10 +479,26 @@ _AssociativeList*    _TreeTopology::MainTreeConstructor  (_String const& parms,
                     }
                     nodeValue = parms.Cut (numerical_match(0), numerical_match(1));
                     i = numerical_match(1);*/
+                    
+                    /** SLKP 20200813
+                        sscanf is going to call strlen every time, which on long strings could be QUITE expensive!
+                        to limit the cost of this, we are going to assume that the float, if it's there is present within kBLLookahead characters of the current position, copy out the next kBLLookahead chars, and read from there.
+                    */
+                    
                     int end_at;
                     hyFloat length = 0.;
-                    if (sscanf (parms.get_str()+i+1, "%lf%n", &length, &end_at) != 1) {
-                        throw _String("Failed to read a number for the branch length following ':'");
+
+                    if (i + kBLLookahead < parms.length()) {
+                        char buffer [kBLLookahead+1];
+                        memcpy (buffer, parms.get_str() + i + 1, kBLLookahead);
+                        buffer [kBLLookahead] = 1;
+                        if (sscanf (buffer, "%lf%n", &length, &end_at) != 1) {
+                            throw _String("Failed to read a number for the branch length following ':'");
+                        }
+                    } else {
+                        if (sscanf (parms.get_str()+i+1, "%lf%n", &length, &end_at) != 1) {
+                            throw _String("Failed to read a number for the branch length following ':'");
+                        }
                     }
                     nodeValue = parms.Cut (i+1, i+end_at);
                     i += end_at;
diff --git a/src/core/tree.cpp b/src/core/tree.cpp
index abf11b939..b5cf314b8 100644
--- a/src/core/tree.cpp
+++ b/src/core/tree.cpp
@@ -2624,42 +2624,36 @@ long    _TheTree::ComputeReleafingCostChar (_DataSetFilter const* dsf, long firs
     const char *pastState = dsf->GetColumn(firstIndex),
                *thisState = dsf->GetColumn(secondIndex);
 
+    long rootIndex = flatTree.lLength - 1,
+         theCost = child_count->list_data[rootIndex],
+         offset = flatLeaves.countitems();
 
     bool * marked_nodes = (bool*) alloca (sizeof(bool) * flatTree.lLength);
     InitializeArray(marked_nodes, flatTree.lLength, false);
-    
-    flatLeaves.Each ([&] (long node, unsigned long node_index) -> void {
-        long f = dsf->theNodeMap.list_data[node_index];
-        if (thisState[f] != pastState[f]) {
-            marked_nodes [flatParents.get(node_index)] = true;
-        }
-    });
-
-    long theCost = 0L,
-         offset = flatLeaves.countitems();
-    
+     
+    for (long node_index = 0L; node_index < flatLeaves.lLength; node_index++) {
+         long seq_index = dsf->theNodeMap.list_data[node_index];
+         if (thisState[seq_index] != pastState[seq_index]) {
+             /*long parentIndex = flatParents.list_data[node_index];
+             while (marked_nodes[parentIndex] == false && parentIndex < rootIndex) {
+                 marked_nodes[parentIndex]  = true;
+                 theCost += child_count->list_data [parentIndex];
+                 parentIndex = flatParents.list_data[parentIndex+offset];
+             }*/
+             marked_nodes [flatParents.list_data[node_index]] = true;
+         }
+    }
 
-    for (long i=0; i<flatTree.lLength; i++) {
+    // the root will always be changed, so don't need to check it
+    for (long i=0; i<flatTree.lLength - 1; i++) {
         if (marked_nodes[i]) {
-            long myParent = flatParents.get (i + offset);
-            if (myParent >= 0) {
-                marked_nodes [myParent] = true;
-            }
-            theCost += child_count->get (i);
+            marked_nodes [flatParents.list_data [i + offset]] = true;
+            theCost += child_count->list_data [i];
         }
     }
 
     return theCost;
-    
-    /*_SimpleList store_nodes (flatTree.lLength, (long*)alloca (flatTree.lLength * sizeof (long)));
-    _AVLList    touched_nodes (&store_nodes);
-    
-    flatLeaves.Each ([&] (long node, unsigned long node_index) -> void {
-        long f = dsf->theNodeMap.list_data[node_index];
-        if (thisState[f] != pastState[f]) {
-            touched_nodes.InsertNumber(flatParents.get(node_index));
-        }
-    });*/
+
 }
 
 //_______________________________________________________________________________________________
@@ -2674,7 +2668,7 @@ void    _TheTree::ClearConstraints (void) {
 
 //_______________________________________________________________________________________________
 
-long    _TheTree::ComputeReleafingCost (_DataSetFilter const* dsf, long firstIndex, long secondIndex, _SimpleList* traversalTags, long orderIndex) const {
+long    _TheTree::ComputeReleafingCost (_DataSetFilter const* dsf, long firstIndex, long secondIndex, _SimpleList* traversalTags, long orderIndex, _SimpleList const* childCount) const {
 
     long        filterL = dsf->GetPatternCount();
 
@@ -2688,22 +2682,44 @@ long    _TheTree::ComputeReleafingCost (_DataSetFilter const* dsf, long firstInd
         marked_nodes [this->flatParents.list_data[di]] = true;
     });
     
-    long theCost = 0;
-
-
-    for (long i=0; i<flatTree.lLength; i++) {
-        if (marked_nodes[i]) {
-            long myParent    =  flatParents.list_data[flatLeaves.lLength + i];
-            if (myParent >= 0) {
-                marked_nodes[myParent] = true;
+    long theCost = 0,
+         rootIndex = flatTree.lLength - 1;
+
+    // don't check the root; it is always "dirty"
+    
+    if (childCount) {
+        if (traversalTags && orderIndex) {
+            for (long i=0; i<rootIndex; i++) {
+               if (marked_nodes[i]) {
+                   marked_nodes[flatParents.list_data[flatLeaves.lLength + i]] = true;
+                   theCost += childCount->list_data[i];
+               } else {
+                   long theIndex = filterL * i + orderIndex;
+                   traversalTags->list_data[theIndex/_HY_BITMASK_WIDTH_] |= bitMaskArray.masks[theIndex%_HY_BITMASK_WIDTH_];
+               }
+           }
+        }
+        
+        for (long i=0; i<rootIndex; i++) {
+            if (marked_nodes[i]) {
+                marked_nodes[flatParents.list_data[flatLeaves.lLength + i]] = true;
+                theCost += childCount->list_data[i];
+            }
+        }
+    } else {
+        for (long i=0; i<rootIndex; i++) {
+            if (marked_nodes[i]) {
+                marked_nodes[flatParents.list_data[flatLeaves.lLength + i]] = true;
+                theCost += ((node <long>*)(flatNodes.list_data[i]))->get_num_nodes();
+            } else if (traversalTags && orderIndex) {
+                long theIndex = filterL * i + orderIndex;
+                traversalTags->list_data[theIndex/_HY_BITMASK_WIDTH_] |= bitMaskArray.masks[theIndex%_HY_BITMASK_WIDTH_];
             }
-            theCost         += ((node <long>*)(flatNodes.list_data[i]))->get_num_nodes();
-        } else if (traversalTags && orderIndex) {
-            long theIndex = filterL * i + orderIndex;
-            traversalTags->list_data[theIndex/_HY_BITMASK_WIDTH_] |= bitMaskArray.masks[theIndex%_HY_BITMASK_WIDTH_];
         }
     }
 
+    theCost += ((node <long>*)(flatNodes.list_data[rootIndex]))->get_num_nodes();
+    
     return theCost;
 
 }
@@ -3154,8 +3170,7 @@ hyFloat      _TheTree::ComputeTreeBlockByBranch  (                   _SimpleList
             
             if (alphabetDimension == 4UL) {
                 long k3     = 0;
-                if (matchSet)
-                    
+                if (matchSet) {
                     for (long k = siteFrom; k < siteTo; k++, k3+=4) {
                         parentConditionals [k3]   = 0.;
                         parentConditionals [k3+1] = 0.;
@@ -3163,14 +3178,18 @@ hyFloat      _TheTree::ComputeTreeBlockByBranch  (                   _SimpleList
                         parentConditionals [k3+3] = 0.;
                         parentConditionals [k3+setBranchTo[siteOrdering.list_data[k]]] = localScalingFactor[k];
                     }
-                else {
-                    
-                    for (long k = siteFrom; k < siteTo; k++, k3+=4) {
+                } else {
+                     for (long k = siteFrom; k < siteTo; k++, k3+=4) {
                         hyFloat scaler = localScalingFactor[k];
                         parentConditionals [k3]   = scaler;
                         parentConditionals [k3+1] = scaler;
                         parentConditionals [k3+2] = scaler;
                         parentConditionals [k3+3] = scaler;
+                         
+                        /*if (parentCode == 32194L && k == 1205L) {
+                        #pragma omp critical
+                            printf ("Site %ld, scaler %lg\n", k, scaler);
+                        }*/
                     }
                 }
             } else {
@@ -3181,8 +3200,43 @@ hyFloat      _TheTree::ComputeTreeBlockByBranch  (                   _SimpleList
                          pp[setBranchTo[siteOrdering.list_data[k]]] = localScalingFactor[k];
                     }
                 } else {
-                    for (long k = siteFrom; k < siteTo; k++, pp += alphabetDimension) {
-                        InitializeArray(pp, alphabetDimension, (hyFloat)localScalingFactor[k]);
+                    if (alphabetDimensionmod4 == 60) {
+                        for (long k = siteFrom; k < siteTo; k++, pp += alphabetDimension) {
+                            hyFloat lsf = localScalingFactor[k];
+                            pp[0] = lsf;pp[1] = lsf;pp[2] = lsf;pp[3] = lsf;
+                            pp[4] = lsf;pp[5] = lsf;pp[6] = lsf;pp[7] = lsf;
+                            pp[8] = lsf;pp[9] = lsf;pp[10] = lsf;pp[11] = lsf;
+                            pp[12] = lsf;pp[13] = lsf;pp[14] = lsf;pp[15] = lsf;
+                            pp[16] = lsf;pp[17] = lsf;pp[18] = lsf;pp[19] = lsf;
+                            pp[20] = lsf;pp[21] = lsf;pp[22] = lsf;pp[23] = lsf;
+                            pp[24] = lsf;pp[25] = lsf;pp[26] = lsf;pp[27] = lsf;
+                            pp[28] = lsf;pp[29] = lsf;pp[30] = lsf;pp[31] = lsf;
+                            pp[32] = lsf;pp[33] = lsf;pp[34] = lsf;pp[35] = lsf;
+                            pp[36] = lsf;pp[37] = lsf;pp[38] = lsf;pp[39] = lsf;
+                            pp[40] = lsf;pp[41] = lsf;pp[42] = lsf;pp[43] = lsf;
+                            pp[44] = lsf;pp[45] = lsf;pp[46] = lsf;pp[47] = lsf;
+                            pp[48] = lsf;pp[49] = lsf;pp[50] = lsf;pp[51] = lsf;
+                            pp[52] = lsf;pp[53] = lsf;pp[54] = lsf;pp[55] = lsf;
+                            pp[56] = lsf;pp[57] = lsf;pp[58] = lsf;pp[59] = lsf;
+                            for (long k2 = alphabetDimensionmod4; k2 < alphabetDimension; k2++) {
+                                pp[k2] = lsf;
+                            }
+                        }
+                    } else {
+                        if (alphabetDimension == 20UL) {
+                            for (long k = siteFrom; k < siteTo; k++, pp += alphabetDimension) {
+                                hyFloat lsf = localScalingFactor[k];
+                                pp[0] = lsf; pp[1] = lsf; pp[2] = lsf; pp[3] = lsf;
+                                pp[4] = lsf; pp[5] = lsf; pp[6] = lsf; pp[7] = lsf;
+                                pp[8] = lsf; pp[9] = lsf; pp[10] = lsf; pp[11] = lsf;
+                                pp[12] = lsf; pp[13] = lsf; pp[14] = lsf; pp[15] = lsf;
+                                pp[16] = lsf; pp[17] = lsf; pp[18] = lsf; pp[19] = lsf;
+                            }
+                        } else {
+                            for (long k = siteFrom; k < siteTo; k++, pp += alphabetDimension) {
+                                InitializeArray(pp, alphabetDimension, (hyFloat)localScalingFactor[k]);
+                            }
+                        }
                     }
                 }
             }
@@ -3256,7 +3310,7 @@ hyFloat      _TheTree::ComputeTreeBlockByBranch  (                   _SimpleList
             hyFloat  const *tMatrix = transitionMatrix;
             hyFloat  sum         = 0.0;
             
-            char        didScale = 0;
+            long        didScale = 0;
             
             if (isLeaf) {
                 long siteState;
@@ -3335,6 +3389,8 @@ hyFloat      _TheTree::ComputeTreeBlockByBranch  (                   _SimpleList
 
             if (alphabetDimension == 4L) { // special case for nuc data
                 
+                //hyFloat stash [4] = {parentConditionals[0], parentConditionals[1], parentConditionals[2], parentConditionals[3]};
+                
 #ifdef _SLKP_USE_AVX_INTRINSICS
                 _handle4x4_pruning_case (childVector, tMatrix, parentConditionals, tmatrix_transpose);
 #else
@@ -3345,8 +3401,29 @@ hyFloat      _TheTree::ComputeTreeBlockByBranch  (                   _SimpleList
                 // for example if a change must happen on a zero-branch length
                 
                 sum     = (parentConditionals [0] + parentConditionals [1]) + (parentConditionals [2] + parentConditionals [3]);
+                
+                /*if (isnan(sum)) {
+                    char buffer[512];
+                    long child_count = 0;
+                    for (long k = 0; k < flatParents.lLength; k++) {
+                        child_count += flatParents.list_data [k] == parentCode;
+                    }
+                    
+                    snprintf (buffer, 512, "Encountered a NAN in parentConditionals at node %s (%ld code, %ld children of parent), site %ld: %lg %lg %lg %lg (%lg, %lg, %lg, %lg); scaler = %lg", currentTreeNode->GetName()->get_str(), parentCode, child_count, siteID, parentConditionals[0],parentConditionals[1], parentConditionals[2], parentConditionals[3],
+                              childVector[0],childVector[1],childVector[2],childVector[3],_lfScalerUpwards);
+                    HandleApplicationError(buffer);
+                }*/
+                
                 if (sum < _lfScalingFactorThreshold && sum > 0.0) {
                     hyFloat tryScale                                 = scalingAdjustments [parentCode*siteCount + siteID] * _lfScalerUpwards;
+                    
+                    /*if (siteID == 1205) {
+                        #pragma omp critical
+                        {
+                            printf ("Scaling UP at node %s, Sum = %lg, scaler = %lg\n",  currentTreeNode->GetName()->get_str(), sum, tryScale);
+                        }
+                    }*/
+                    
                     if (tryScale < HUGE_VAL) {
                         parentConditionals [0]                             *= _lfScalerUpwards;
                         parentConditionals [1]                             *= _lfScalerUpwards;
@@ -3358,14 +3435,36 @@ hyFloat      _TheTree::ComputeTreeBlockByBranch  (                   _SimpleList
                     }
                 } else {
                     if (sum > _lfScalerUpwards) {
-                        scalingAdjustments [parentCode*siteCount + siteID] *= _lfScalingFactorThreshold;
-                        parentConditionals [0]                             *= _lfScalingFactorThreshold;
-                        parentConditionals [1]                             *= _lfScalingFactorThreshold;
-                        parentConditionals [2]                             *= _lfScalingFactorThreshold;
-                        parentConditionals [3]                             *= _lfScalingFactorThreshold;
-                        localScalerChange                                  -= theFilter->theFrequencies.get (siteOrdering.list_data[siteID]);
-                        didScale                                            = -1;
+                        if (sum < HUGE_VAL) { // no point scaling an infinity
+                            
+                            didScale                                            = -1;
+                            
+                            hyFloat scaler                                 = scalingAdjustments [parentCode*siteCount + siteID] * _lfScalingFactorThreshold;
+                            //hyFloat save_sum = sum;
+                            sum *= _lfScalingFactorThreshold;
+                            while (sum > _lfScalerUpwards) {
+                                sum     *= _lfScalingFactorThreshold;
+                                scaler  *= _lfScalingFactorThreshold;
+                                didScale --;
+                            }
+                            
+                            /*#pragma omp critical
+                            if (siteID == 1205) {
+                                printf ("Scaling down from %lg to %lg; scaler %lg, steps %ld\n", save_sum, sum, scaler, didScale);
+                            }*/
+                            
+                            parentConditionals [0]                             *= scaler;
+                            parentConditionals [1]                             *= scaler;
+                            parentConditionals [2]                             *= scaler;
+                            parentConditionals [3]                             *= scaler;
+                           
+                            localScalerChange                                  += didScale * theFilter->theFrequencies (siteOrdering.list_data[siteID]);
+                            scalingAdjustments [parentCode*siteCount + siteID] = scaler;
+                            //didScale                                         = -1;
+                        }
+                        
                     }
+                
                 }
                 
                 childVector += 4;
@@ -3487,7 +3586,7 @@ hyFloat      _TheTree::ComputeTreeBlockByBranch  (                   _SimpleList
 #endif // regular code
                         
                         if (alphabetDimension == 61) {
-                            sum += (parentConditionals[p] *= accumulator + tMatrix[alphabetDimensionmod4] * childVector[alphabetDimensionmod4]);
+                            sum += (parentConditionals[p] *= accumulator + tMatrix[60] * childVector[60]);
                         } else {
                             for (long c = alphabetDimensionmod4; c < alphabetDimension; c++) {
                                 accumulator +=  tMatrix[c] * childVector[c];
@@ -3635,20 +3734,23 @@ hyFloat      _TheTree::ComputeTreeBlockByBranch  (                   _SimpleList
             long                rootState = setBranchTo[siteOrdering.list_data[siteID]];
             accumulator         = rootConditionals[rootIndex + rootState] * theProbs[rootState];
             rootIndex           += alphabetDimension;
-        } else
+        } else {
+            
             for (long p = 0; p < alphabetDimension; p++,rootIndex++) {
                 accumulator += rootConditionals[rootIndex] * theProbs[p];
+                
             }
+        }
+           
+        /*#pragma omp critical
+        if (siteID == 1205) {
+            printf ("\n ACCUMULATOR = %lg (%ld)\n", accumulator, likeFuncEvalCallCount);
+        }*/
                 
         /*#pragma omp critical
          {
-         if (likeFuncEvalCallCount == 0) {
-         eval_buffer [siteID] = accumulator;
-         } else {
-         if (!CheckEqual (eval_buffer [siteID], accumulator)) {
-         printf ("EVAL %ld, Site %ld/%ld, LogL %g, %g (%d)\n", likeFuncEvalCallCount, siteFrom, siteID, accumulator, eval_buffer [siteID],localScalerChange);
-         }
-         eval_buffer [siteID] = accumulator;
+         if (true || likeFuncEvalCallCount == 8) {
+             printf ("EVAL %ld, Site %ld/%ld, %g (freq = %d, log L = %g), (%g) is le0 = %d\n", likeFuncEvalCallCount, siteID, siteTo, accumulator , theFilter->theFrequencies (siteOrdering.list_data[siteID]), log (accumulator), correction, accumulator <= 0.0);
          }
          }*/
         
@@ -3659,6 +3761,7 @@ hyFloat      _TheTree::ComputeTreeBlockByBranch  (                   _SimpleList
                 result = -INFINITY;
 #pragma omp critical
                 {
+                    //printf ("BAILING WITH INFINITY %ld\n", siteID);
                     hy_global::ReportWarning (_String("Site ") & (1L+siteOrdering.list_data[siteID]) & " evaluated to a 0 probability in ComputeTreeBlockByBranch");
                 }
                 break;
@@ -3965,7 +4068,7 @@ void            _TheTree::ComputeBranchCache    (
             }
             
             hyFloat sum      = .0;
-            char       didScale = 0;
+            long       didScale = 0;
             
             if (alphabetDimension == 4UL) { // special case for nuc data
 #ifdef _SLKP_USE_AVX_INTRINSICS
@@ -3986,16 +4089,35 @@ void            _TheTree::ComputeBranchCache    (
                             
                             localScalerChange                                  += theFilter->theFrequencies (siteOrdering.list_data[siteID]);
                             didScale                                            = 1;
+                            scalingAdjustments [nodeCode*siteCount + siteID]    = tryScale;
                         }
                     } else {
                         if (sum > _lfScalerUpwards) {
-                            parentConditionals [0]                             *= _lfScalingFactorThreshold;
-                            parentConditionals [1]                             *= _lfScalingFactorThreshold;
-                            parentConditionals [2]                             *= _lfScalingFactorThreshold;
-                            parentConditionals [3]                             *= _lfScalingFactorThreshold;
+                            if (sum < HUGE_VAL) { // no point scaling an infinity
+                                
+                                didScale                                            = -1;
+                                
+                                
+                                hyFloat scaler = scalingAdjustments [parentCode*siteCount + siteID] * _lfScalingFactorThreshold;
+                                
+                                sum *= _lfScalingFactorThreshold;
+                                while (sum > _lfScalerUpwards) {
+                                    sum     *= _lfScalingFactorThreshold;
+                                    scaler  *= _lfScalingFactorThreshold;
+                                    didScale --;
+                                }
+                                
+                                
+                                parentConditionals [0]                             *= scaler;
+                                parentConditionals [1]                             *= scaler;
+                                parentConditionals [2]                             *= scaler;
+                                parentConditionals [3]                             *= scaler;
+                               
+                                localScalerChange                                  += didScale * theFilter->theFrequencies (siteOrdering.list_data[siteID]);
+                                //didScale                                         = -1;
+                                scalingAdjustments [nodeCode*siteCount + siteID]    = scaler;
+                            }
                             
-                            localScalerChange                                  -= theFilter->theFrequencies (siteOrdering.list_data[siteID]);
-                            didScale                                            = -1;
                         }
                     }
                 }
diff --git a/src/core/variablecontainer.cpp b/src/core/variablecontainer.cpp
index db476fc8a..ebfcea8d4 100644
--- a/src/core/variablecontainer.cpp
+++ b/src/core/variablecontainer.cpp
@@ -853,6 +853,8 @@ void      _VariableContainer::ClearConstraints(void) {
     }
 }
 
+//#define _UBER_VERBOSE_MX_UPDATE_DUMP
+
 //__________________________________________________________________________________
 
 void      _VariableContainer::CopyModelParameterValue (long var_idx, long ref_idx, unsigned long) {
@@ -862,9 +864,7 @@ void      _VariableContainer::CopyModelParameterValue (long var_idx, long ref_id
             model_var->SetValue (LocateVar (var_idx)->Compute(),true,true,NULL);
             
 #ifdef _UBER_VERBOSE_MX_UPDATE_DUMP
-            if (1 || l == _UBER_VERBOSE_MX_UPDATE_DUMP_LF_EVAL) {
                 fprintf (stderr, "[_CalcNode::RecomputeMatrix] Node %s, var %s, value = %15.12g\n", LocateVar (var_idx)->GetName()->get_str(), model_var->GetName()->get_str(), model_var->Compute()->Value());
-            }
 #endif
         }
     }
diff --git a/src/mains/unix.cpp b/src/mains/unix.cpp
index c675ad894..fa19d1e09 100644
--- a/src/mains/unix.cpp
+++ b/src/mains/unix.cpp
@@ -691,7 +691,8 @@ void    SetStatusLineUser   (_String const s) {
 
 #ifndef __UNITTEST__
 int main (int argc, char* argv[]) {
-
+    
+ 
 #ifdef _COMPARATIVE_LF_DEBUG_DUMP
     FILE * comparative_lf_debug_matrix_content_file = doFileOpen (_COMPARATIVE_LF_DEBUG_DUMP, "w");
 #endif
diff --git a/tests/hbltests/Results/HIVSweden.out b/tests/hbltests/Results/HIVSweden.out
index 127c7e7b2..79d8209f6 100644
--- a/tests/hbltests/Results/HIVSweden.out
+++ b/tests/hbltests/Results/HIVSweden.out
@@ -4,11 +4,11 @@
 
 >Done in 4 seconds 
 
--1137.68851199504
+-1137.68851199166
 
 -----------------------------------
 
-dN/dS = 0.9053687405962891
+dN/dS = 0.9053687083867886
 
 
 *** RUNNING MODEL 1 (Neutral) ***
@@ -16,14 +16,14 @@ dN/dS = 0.9053687405962891
 
 >Done in 2 seconds 
 
--1114.64186455482
+-1114.64201891122
 
 ------------------------------------------------
 
-dN/dS = 0.554584342157782 (sample variance = 0.2120936670031687)
+dN/dS = 0.5549352383990346 (sample variance = 0.2120727849648286)
 
-Rate[1]=  0.07841413 (weight=0.4833143)
-Rate[2]=  1.00000000 (weight=0.5166857)
+Rate[1]=  0.07843653 (weight=0.4829453)
+Rate[2]=  1.00000000 (weight=0.5170547)
 
 ------------------------------------------------
 
@@ -37,215 +37,215 @@ Import the following part into a data processing program
 for further analysis
 
 
-Rate/Site	Rate=0.07841413289607745	Rate=1
-1	0.000413882148180681	0.9995861178518194
-2	0.9496498025372991	0.05035019746270093
-3	0.8530116627311244	0.1469883372688757
-4	0.9773824962278876	0.02261750377211239
-5	0.7092729013090043	0.2907270986909957
-6	0.9707183515789	0.02928164842109998
-7	0.8414170785332546	0.1585829214667454
-8	0.6451742850403653	0.3548257149596347
-9	0.009009098906633749	0.9909909010933663
-10	0.2516636715151961	0.7483363284848039
-11	0.9364433956371618	0.06355660436283812
-12	0.8043655892037351	0.1956344107962648
-13	0.9594118973469086	0.04058810265309146
-14	0.6835043668332426	0.3164956331667573
-15	0.968787475058735	0.03121252494126506
-16	0.968787475058735	0.03121252494126506
-17	0.9496498025372991	0.05035019746270093
-18	0.0962431376466097	0.9037568623533904
-19	0.8049051987122308	0.1950948012877692
-20	0.4256863428153749	0.574313657184625
-21	0.00115841682150587	0.9988415831784941
-22	0.001176884510908392	0.9988231154890916
-23	0.9496498025372991	0.05035019746270093
-24	7.029769746973078e-05	0.9999297023025302
-25	0.8603172169344939	0.1396827830655061
-26	2.227461563248114e-05	0.9999777253843676
-27	0.716145078807327	0.2838549211926729
-28	1.803396035103997e-08	0.9999999819660396
-29	0.9786676818541649	0.02133231814583518
-30	0.5827249653652521	0.4172750346347479
-31	0.002372181619193272	0.9976278183808068
-32	0.8414170785332546	0.1585829214667454
-33	0.9364433956371618	0.06355660436283812
-34	0.698370682494158	0.301629317505842
-35	0.7788091406650403	0.2211908593349597
-36	0.01418626148588906	0.9858137385141109
-37	0.06334915401471131	0.9366508459852887
-38	0.6916538615990294	0.3083461384009705
-39	0.001913788952689397	0.9980862110473107
-40	0.003952551303134973	0.9960474486968649
-41	0.5299501283259098	0.4700498716740902
-42	0.8386599846750767	0.1613400153249232
-43	0.5925549851664365	0.4074450148335635
-44	0.2083179670582381	0.791682032941762
-45	0.4815786235964578	0.5184213764035422
-46	0.01852084337622002	0.9814791566237799
-47	0.7641444205947111	0.2358555794052888
-48	0.140302665422494	0.859697334577506
-49	0.07638728509206009	0.9236127149079398
-50	0.5570428556565402	0.4429571443434598
-51	4.373758390987692e-06	0.999995626241609
-52	0.968787475058735	0.03121252494126506
-53	0.9340205870835026	0.06597941291649743
-54	0.05710779684898481	0.9428922031510152
-55	0.6737632882461265	0.3262367117538734
-56	0.9786676818541649	0.02133231814583518
-57	0.03633265439224938	0.9636673456077506
-58	0.9368002800854486	0.06319971991455135
-59	0.02359697208779809	0.9764030279122018
-60	0.716145078807327	0.2838549211926729
-61	0.2740753812362675	0.7259246187637325
-62	0.1720713170353596	0.8279286829646405
-63	0.06693254813189572	0.9330674518681042
-64	0.1591604900838105	0.8408395099161895
-65	0.1371356664179075	0.8628643335820925
-66	2.815619623961706e-08	0.9999999718438037
-67	0.5393523345882355	0.4606476654117646
-68	0.0006606580284920114	0.999339341971508
-69	2.85469434309379e-05	0.9999714530565691
-70	0.3236625620665379	0.6763374379334621
-71	0.2249319499347314	0.7750680500652686
-72	0.1619254042736656	0.8380745957263342
-73	0.1394106370415928	0.860589362958407
-74	0.9496498025372991	0.05035019746270093
-75	0.02456229364895695	0.975437706351043
-76	6.382038068508085e-05	0.999936179619315
-77	0.9716280237009275	0.02837197629907249
-78	0.7865549915690062	0.2134450084309938
-79	0.7333992239855914	0.2666007760144085
-80	0.9707183515789	0.02928164842109998
-81	0.8613584674737541	0.1386415325262458
-82	0.7671225351999779	0.2328774648000221
-83	5.644658708439183e-05	0.9999435534129156
-84	0.006716287357999313	0.9932837126420007
-85	0.9684714958300438	0.03152850416995622
-86	0.968787475058735	0.03121252494126506
-87	1.515713382192203e-05	0.9999848428661781
-88	0.7671225351999779	0.2328774648000221
-89	0.9317677181261299	0.06823228187387004
-90	0.03855611577601334	0.9614438842239867
-91	0.5875486193836759	0.4124513806163241
+Rate/Site	Rate=0.07843652795974788	Rate=1
+1	0.0004129669814584291	0.9995870330185416
+2	0.9494095652569839	0.05059043474301613
+3	0.8526035606507344	0.1473964393492656
+4	0.9772565198354999	0.02274348016450016
+5	0.7087876841962609	0.291212315803739
+6	0.9705564425354483	0.02944355746455165
+7	0.8410477831543381	0.1589522168456618
+8	0.6447244792865031	0.3552755207134968
+9	0.008996956992291881	0.9910030430077081
+10	0.2511779934992873	0.7488220065007126
+11	0.9362195651764531	0.06378043482354691
+12	0.8038489724103089	0.1961510275896911
+13	0.959236190913887	0.04076380908611305
+14	0.6829847717701353	0.3170152282298647
+15	0.9686391929170401	0.03136080708295985
+16	0.9686391929170401	0.03136080708295985
+17	0.9494095652569839	0.05059043474301613
+18	0.09612255426779577	0.9038774457322042
+19	0.8045219801863088	0.1954780198136912
+20	0.424877156861682	0.5751228431383181
+21	0.001154443467410566	0.9988455565325894
+22	0.00117506086167681	0.9988249391383233
+23	0.9494095652569839	0.05059043474301613
+24	7.003513198056934e-05	0.9999299648680195
+25	0.8599134733983664	0.1400865266016336
+26	2.223230551768934e-05	0.9999777676944822
+27	0.7156452019595485	0.2843547980404515
+28	1.803036482298124e-08	0.9999999819696351
+29	0.9785466252430841	0.02145337475691597
+30	0.5822748962818075	0.4177251037181926
+31	0.002367566727977607	0.9976324332720224
+32	0.8410477831543381	0.1589522168456618
+33	0.9362195651764531	0.06378043482354691
+34	0.6973646326686107	0.3026353673313892
+35	0.7778580020052748	0.2221419979947252
+36	0.01413624080003053	0.9858637591999695
+37	0.06314701285732313	0.936852987142677
+38	0.6906729097755763	0.3093270902244237
+39	0.001910666259187429	0.9980893337408127
+40	0.003938810644126677	0.9960611893558734
+41	0.5287725420386381	0.471227457961362
+42	0.8382470749672665	0.1617529250327335
+43	0.5920941790931827	0.4079058209068173
+44	0.2074662887277795	0.7925337112722205
+45	0.4806045408976926	0.5193954591023074
+46	0.01848817363411778	0.9815118263658823
+47	0.7636899411185869	0.236310058881413
+48	0.139769181967181	0.8602308180328189
+49	0.07626105691443552	0.9237389430855645
+50	0.5558696207398601	0.4441303792601399
+51	4.361524767754885e-06	0.9999956384752322
+52	0.9686391929170401	0.03136080708295985
+53	0.9337130467062972	0.06628695329370292
+54	0.05692188471752498	0.943078115282475
+55	0.6727699939355163	0.3272300060644837
+56	0.9785466252430841	0.02145337475691597
+57	0.03625726463428663	0.9637427353657133
+58	0.9365169892139765	0.06348301078602352
+59	0.02354995388334547	0.9764500461166545
+60	0.7156452019595485	0.2843547980404515
+61	0.2736138415085682	0.7263861584914317
+62	0.1717726094741855	0.8282273905258145
+63	0.06671376548001214	0.9332862345199878
+64	0.1584764771790627	0.8415235228209373
+65	0.1366678299505407	0.8633321700494593
+66	2.806195822890807e-08	0.9999999719380418
+67	0.5389526198273991	0.461047380172601
+68	0.0006594279056825333	0.9993405720943175
+69	2.851527490764395e-05	0.9999714847250923
+70	0.3229584874662094	0.6770415125337907
+71	0.224552996956325	0.775447003043675
+72	0.1613487030553297	0.8386512969446702
+73	0.1392125323982454	0.8607874676017546
+74	0.9494095652569839	0.05059043474301613
+75	0.02449671329197509	0.9755032867080249
+76	6.360145714514582e-05	0.9999363985428549
+77	0.9715055205623275	0.0284944794376725
+78	0.7861409694662842	0.2138590305337158
+79	0.7323270300560397	0.2676729699439603
+80	0.9705564425354483	0.02944355746455165
+81	0.8610502888731348	0.1389497111268651
+82	0.7666818522524143	0.2333181477475857
+83	5.6313393612785e-05	0.9999436866063872
+84	0.006701121812101592	0.9932988781878984
+85	0.9683170085554258	0.03168299144457424
+86	0.9686391929170401	0.03136080708295985
+87	1.524689710895521e-05	0.9999847531028911
+88	0.7666818522524143	0.2333181477475857
+89	0.9315263734615309	0.06847362653846911
+90	0.03847784219918648	0.9615221578008135
+91	0.5870967803100652	0.4129032196899347
 
 *** RUNNING MODEL 2 (Selection) ***
 ######################################
 
->Done in 8 seconds 
+>Done in 7 seconds 
 
--1106.44526776485
+-1106.4453042489
 
 ------------------------------------------------
 
-dN/dS = 1.12655619184039 (sample variance = 1.580581693007084)
+dN/dS = 1.126982070134336 (sample variance = 1.581623693821516)
 
-Rate[1]=  0.05881756 (weight=0.3752209)
-Rate[2]=  1.00000000 (weight=0.4427542)
-Rate[3]=  3.63539426 (weight=0.1820249)
+Rate[1]=  0.05836173 (weight=0.3746328)
+Rate[2]=  1.00000000 (weight=0.4435036)
+Rate[3]=  3.63796971 (weight=0.1818636)
 
 ------------------------------------------------
 
  Sites with dN/dS>1 (Posterior cutoff = 0.9)
 
-26 (0.9064573245737738)
-28 (0.9992020791009154)
-66 (0.9984629071680293)
-87 (0.9855170546898564)
+26 (0.9063116682671567)
+28 (0.9992009807387726)
+66 (0.9984632592060946)
+87 (0.9855127564740064)
 
 
 ------------------------------------------------
 
  Sites with dN/dS<=1 (Posterior cutoff = 0.9)
 
-1 (0.5037453402632631)
-2 (0.0004588106769083485)
-3 (0.005906139478386191)
-4 (9.654170642798585e-05)
-5 (0.03166315408305304)
-6 (0.000151012915826745)
-7 (0.006934491590164102)
-8 (0.05661370569906062)
-9 (0.7216327366457385)
-10 (0.09628729296561687)
-11 (0.0009199027526116262)
-12 (0.01061995625253789)
-13 (0.0004927618917987312)
-14 (0.03982762843583483)
-15 (0.0002223182636933575)
-16 (0.0002223182636933575)
-17 (0.0004588106769083485)
-18 (0.4321134585058325)
-19 (0.01159547152948417)
-20 (0.03360117046826709)
-21 (0.2465600025243685)
-22 (0.7975622136581251)
-23 (0.0004588106769083485)
-24 (0.5827716819681953)
-25 (0.005164666146267199)
-27 (0.02934394809734361)
-29 (7.94606599295118e-05)
-30 (0.09124445227479458)
-31 (0.5688123449123863)
-32 (0.006934491590164102)
-33 (0.0009199027526116262)
-34 (0.004506156393895435)
-35 (0.002046261317851369)
-36 (0.104959954770928)
-37 (0.08582024673150629)
-38 (0.005386803621285784)
-39 (0.6429273833262951)
-40 (0.4544859505309178)
-41 (0.01182754218631978)
-42 (0.006847995833739094)
-43 (0.08347425154688627)
-44 (0.01182353194578844)
-45 (0.01747472136380717)
-46 (0.4284692663941158)
-47 (0.0183451423381276)
-48 (0.02531160529659541)
-49 (0.07103029501470552)
-50 (0.0106118372575389)
-51 (0.8853156372325378)
-52 (0.0002223182636933575)
-53 (0.0007673192996036272)
-54 (0.09813360180044658)
-55 (0.006327000138082545)
-56 (7.94606599295118e-05)
-57 (0.2046375777641244)
-58 (0.0007266494606770322)
-59 (0.3307042574250877)
-60 (0.02934394809734361)
-61 (0.08182417232921647)
-62 (0.1931989415880896)
-63 (0.07709868732161493)
-64 (0.02144654836661221)
-65 (0.02740837592622346)
-67 (0.1258714465550827)
-68 (0.6042881515947164)
-69 (0.8324348912682272)
-70 (0.05521683558757644)
-71 (0.1225064221507603)
-72 (0.0210794443186693)
-73 (0.2770274713029749)
-74 (0.0004588106769083485)
-75 (0.3065817804707144)
-76 (0.6746029184227686)
-77 (0.0002302421435441318)
-78 (0.01496861044592552)
-79 (0.003216410716239359)
-80 (0.000151012915826745)
-81 (0.005168660879972707)
-82 (0.01773982794920307)
-83 (0.8132652038617769)
-84 (0.2778069717762056)
-85 (0.0001963012176694039)
-86 (0.0002223182636933575)
-88 (0.01773982794920307)
-89 (0.001171879566541098)
-90 (0.1860446585703643)
-91 (0.08807631966319565)
+1 (0.5029882717724251)
+2 (0.0004565934951104923)
+3 (0.005887938153177528)
+4 (9.600689888168994e-05)
+5 (0.03159578682248793)
+6 (0.0001501912317134283)
+7 (0.006913595036040512)
+8 (0.05651627369783895)
+9 (0.7213059006955111)
+10 (0.09603380485043612)
+11 (0.0009161086971956588)
+12 (0.0105892160733136)
+13 (0.0004906877706081168)
+14 (0.03974830839056758)
+15 (0.0002212421990611336)
+16 (0.0002212421990611336)
+17 (0.0004565934951104923)
+18 (0.4316844811961461)
+19 (0.01156424576429309)
+20 (0.03349647980970066)
+21 (0.2457975252198305)
+22 (0.7972570065707351)
+23 (0.0004565934951104923)
+24 (0.5819108177701073)
+25 (0.005148229652313554)
+27 (0.02927961337167887)
+29 (7.90015239951085e-05)
+30 (0.09111607918977237)
+31 (0.5680980118337913)
+32 (0.006913595036040512)
+33 (0.0009161086971956588)
+34 (0.004491013525954807)
+35 (0.002039198916549559)
+36 (0.104507972231377)
+37 (0.08547148167400613)
+38 (0.005370721667639105)
+39 (0.6423404619343602)
+40 (0.4536876051126541)
+41 (0.01178307858562074)
+42 (0.006826751819136865)
+43 (0.0833467202789617)
+44 (0.01175537509306583)
+45 (0.01741197798098377)
+46 (0.4278075925123913)
+47 (0.01829998147118346)
+48 (0.02518055052153371)
+49 (0.07070853519209956)
+50 (0.01057248191459843)
+51 (0.8850702194936431)
+52 (0.0002212421990611336)
+53 (0.0007637587433421098)
+54 (0.09774658921522028)
+55 (0.006307263649733691)
+56 (7.90015239951085e-05)
+57 (0.2040552822831782)
+58 (0.0007233015487799649)
+59 (0.3300058913730194)
+60 (0.02927961337167887)
+61 (0.08160697871997806)
+62 (0.1928097635222059)
+63 (0.07677075293749337)
+64 (0.02132547110184272)
+65 (0.02727346154992843)
+67 (0.125725778974423)
+68 (0.6035808825405923)
+69 (0.8321722198290895)
+70 (0.05505312153448871)
+71 (0.1222135186938459)
+72 (0.02097187843855246)
+73 (0.2765883351011806)
+74 (0.0004565934951104923)
+75 (0.3058760148284225)
+76 (0.6738839489534654)
+77 (0.0002291803057879414)
+78 (0.01493054914017372)
+79 (0.00320622974990686)
+80 (0.0001501912317134283)
+81 (0.005152396719899931)
+82 (0.01769599223513166)
+83 (0.8129016851732231)
+84 (0.2768950821445459)
+85 (0.0001952974875795364)
+86 (0.0002212421990611336)
+88 (0.01769599223513166)
+89 (0.001167310596224406)
+90 (0.1854866382325809)
+91 (0.08795026739358358)
 
 
 ------------------------------------------------
@@ -256,95 +256,95 @@ Import the following part into a data processing program
 for further analysis
 
 
-Rate/Site	Rate=0.05881756131775031	Rate=1	Rate=3.635394261991776
-1	2.177215825406107e-05	0.4962328875784827	0.5037453402632631
-2	0.8458520763618315	0.1536891129612601	0.0004588106769083485
-3	0.7176381382272181	0.2764557222943958	0.005906139478386191
-4	0.9026014904626679	0.09730196783090417	9.654170642798585e-05
-5	0.5833471670132869	0.3849896789036601	0.03166315408305304
-6	0.8877906508430571	0.1120583362411163	0.000151012915826745
-7	0.7017126840971444	0.2913528243126914	0.006934491590164102
-8	0.5262607791203352	0.4171255151806042	0.05661370569906062
-9	0.001068345191200685	0.277298918163061	0.7216327366457385
-10	0.1010580530030604	0.8026546540313226	0.09628729296561687
-11	0.819785528318149	0.1792945689292394	0.0009199027526116262
-12	0.6722219983714567	0.3171580453760055	0.01061995625253789
-13	0.8649302785570022	0.134576959551199	0.0004927618917987312
-14	0.5628664910358208	0.3973058805283444	0.03982762843583483
-15	0.882531831009213	0.1172458507270937	0.0002223182636933575
-16	0.882531831009213	0.1172458507270937	0.0002223182636933575
-17	0.8458520763618315	0.1536891129612601	0.0004588106769083485
-18	0.03322193659119953	0.5346646049029681	0.4321134585058325
-19	0.6650010053547115	0.3234035231158043	0.01159547152948417
-20	0.1606693511884449	0.805729478343288	0.03360117046826709
-21	6.644705000860961e-05	0.7533735504256228	0.2465600025243685
-22	6.270630524073773e-05	0.2023750800366342	0.7975622136581251
-23	0.8458520763618315	0.1536891129612601	0.0004588106769083485
-24	1.864437782554034e-06	0.4172264535940222	0.5827716819681953
-25	0.7263383110412817	0.2684970228124511	0.005164666146267199
-26	1.946323040306041e-07	0.09354248079392219	0.9064573245737738
-27	0.5900998726138736	0.3805561792887828	0.02934394809734361
-28	7.61853002115326e-14	0.0007979208990084675	0.9992020791009154
-29	0.9057884233713861	0.0941321159686844	7.94606599295118e-05
-30	0.4719272140936588	0.4368283336315467	0.09124445227479458
-31	0.0002016115910458734	0.4309860434965678	0.5688123449123863
-32	0.7017126840971444	0.2913528243126914	0.006934491590164102
-33	0.819785528318149	0.1792945689292394	0.0009199027526116262
-34	0.2976592721127369	0.6978345714933676	0.004506156393895435
-35	0.3577709678773954	0.6401827708047533	0.002046261317851369
-36	0.001331546948636047	0.8937084982804361	0.104959954770928
-37	0.0119508664598319	0.9022288868086619	0.08582024673150629
-38	0.2956738648719643	0.6989393315067498	0.005386803621285784
-39	0.0001468539713052756	0.3569257627023997	0.6429273833262951
-40	0.0002579005312356573	0.5452561489378467	0.4544859505309178
-41	0.2133006053361861	0.7748718524774941	0.01182754218631978
-42	0.7011273459358308	0.2920246582304302	0.006847995833739094
-43	0.4845473455185078	0.431978402934606	0.08347425154688627
-44	0.03048854192233652	0.957687926131875	0.01182353194578844
-45	0.1904137843976104	0.7921114942385824	0.01747472136380717
-46	0.003402221489319509	0.5681285121165648	0.4284692663941158
-47	0.6305177311186616	0.3511371265432107	0.0183451423381276
-48	0.02233164826189267	0.9523567464415119	0.02531160529659541
-49	0.01371463891192875	0.9152550660733657	0.07103029501470552
-50	0.2253661618261848	0.7640220009162764	0.0106118372575389
-51	2.587671966482711e-08	0.1146843368907424	0.8853156372325378
-52	0.882531831009213	0.1172458507270937	0.0002223182636933575
-53	0.8215681550163814	0.1776645256840149	0.0007673192996036272
-54	0.01113272176027104	0.8907336764392825	0.09813360180044658
-55	0.2848448387970111	0.7088281610649063	0.006327000138082545
-56	0.9057884233713861	0.0941321159686844	7.94606599295118e-05
-57	0.007128053604395912	0.7882343686314797	0.2046375777641244
-58	0.8244544813983278	0.1748188691409952	0.0007266494606770322
-59	0.004636723260068558	0.6646590193148437	0.3307042574250877
-60	0.5900998726138736	0.3805561792887828	0.02934394809734361
-61	0.1084510021628918	0.8097248255078917	0.08182417232921647
-62	0.06904784430953879	0.7377532141023716	0.1931989415880896
-63	0.0125729578338336	0.9103283548445515	0.07709868732161493
-64	0.02419770633061857	0.9543557453027692	0.02144654836661221
-65	0.02177593034695595	0.9508156937268206	0.02740837592622346
-66	2.800079631670575e-12	0.001537092829170624	0.9984629071680293
-67	0.4299732415812869	0.4441553118636304	0.1258714465550827
-68	3.149123939333088e-05	0.3956803571658903	0.6042881515947164
-69	4.376693791193806e-07	0.1675646710623936	0.8324348912682272
-70	0.1296615446146847	0.8151216197977387	0.05521683558757644
-71	0.08967799843921676	0.787815579410023	0.1225064221507603
-72	0.0252922101227979	0.9536283455585328	0.0210794443186693
-73	0.05334085012685787	0.6696316785701673	0.2770274713029749
-74	0.8458520763618315	0.1536891129612601	0.0004588106769083485
-75	0.005042518922891693	0.6883757006063939	0.3065817804707144
-76	1.314752652998064e-06	0.3253957668245783	0.6746029184227686
-77	0.8886524865268747	0.1111172713295813	0.0002302421435441318
-78	0.6489825493661842	0.3360488401878903	0.01496861044592552
-79	0.3339539866278762	0.6628296026558844	0.003216410716239359
-80	0.8877906508430571	0.1120583362411163	0.000151012915826745
-81	0.719750934778427	0.2750804043416004	0.005168660879972707
-82	0.6323963707640347	0.3498638012867624	0.01773982794920307
-83	9.724033336775255e-07	0.1867338237348894	0.8132652038617769
-84	9.471843671886234e-05	0.7220983097870755	0.2778069717762056
-85	0.8810742572408299	0.1187294415415008	0.0001963012176694039
-86	0.882531831009213	0.1172458507270937	0.0002223182636933575
-87	1.552866399054062e-08	0.01448292978147969	0.9855170546898564
-88	0.6323963707640347	0.3498638012867624	0.01773982794920307
-89	0.8121298032348194	0.1866983171986396	0.001171879566541098
-90	0.007570515847640617	0.8063848255819951	0.1860446585703643
-91	0.4762812171135835	0.4356424632232208	0.08807631966319565
+Rate/Site	Rate=0.05836173276024324	Rate=1	Rate=3.637969705855394
+1	2.107020468677557e-05	0.4969906580228882	0.5029882717724251
+2	0.8455383320112253	0.1540050744936643	0.0004565934951104923
+3	0.7170945516342836	0.2770175102125389	0.005887938153177528
+4	0.9024086514723629	0.09749534162875539	9.600689888168994e-05
+5	0.582687080902056	0.3857171322754561	0.03159578682248793
+6	0.8875671919014441	0.1122826168668424	0.0001501912317134283
+7	0.7011418419732194	0.2919445629907402	0.006913595036040512
+8	0.5255929470601607	0.4178907792420004	0.05651627369783895
+9	0.001049799103137772	0.277644300201351	0.7213059006955111
+10	0.1001340864333319	0.803832108716232	0.09603380485043612
+11	0.8194118587079032	0.1796720325949011	0.0009161086971956588
+12	0.6716276382027104	0.317783145723976	0.0105892160733136
+13	0.8646623798514339	0.1348469323779579	0.0004906877706081168
+14	0.5622052224963258	0.3980464691131066	0.03974830839056758
+15	0.8822943389422144	0.1174844188587244	0.0002212421990611336
+16	0.8822943389422144	0.1174844188587244	0.0002212421990611336
+17	0.8455383320112253	0.1540050744936643	0.0004565934951104923
+18	0.03290553531167676	0.5354099834921771	0.4316844811961461
+19	0.6643873412055457	0.3240484130301612	0.01156424576429309
+20	0.159293403339795	0.8072101168505044	0.03349647980970066
+21	6.432728885530523e-05	0.7541381474913142	0.2457975252198305
+22	6.116798835512746e-05	0.2026818254409098	0.7972570065707351
+23	0.8455383320112253	0.1540050744936643	0.0004565934951104923
+24	1.792732886383816e-06	0.4180873894970061	0.5819108177701073
+25	0.7258086423329762	0.2690431280147103	0.005148229652313554
+26	1.869722398852567e-07	0.09368814476060341	0.9063116682671567
+27	0.5894432210638682	0.3812771655644528	0.02927961337167887
+28	7.239879007290159e-14	0.000799019261155045	0.9992009807387726
+29	0.9056028772991068	0.09431812117689807	7.90015239951085e-05
+30	0.4712770302444826	0.437606890565745	0.09111607918977237
+31	0.0001967343611938441	0.4317052538050148	0.5680980118337913
+32	0.7011418419732194	0.2919445629907402	0.006913595036040512
+33	0.8194118587079032	0.1796720325949011	0.0009161086971956588
+34	0.2955604193309906	0.6999485671430545	0.004491013525954807
+35	0.3554895577275244	0.642471243355926	0.002039198916549559
+36	0.001298810131237448	0.8941932176373855	0.104507972231377
+37	0.01174275658884314	0.9027857617371506	0.08547148167400613
+38	0.2935827567785841	0.7010465215537768	0.005370721667639105
+39	0.00014328619115538	0.3575162518744844	0.6423404619343602
+40	0.000251580045425019	0.5460608148419208	0.4536876051126541
+41	0.2115982379551061	0.7766186834592732	0.01178307858562074
+42	0.7005588865400794	0.2926143616407837	0.006826751819136865
+43	0.4839037627222449	0.4327495169987934	0.0833467202789617
+44	0.02997078656054276	0.9582738383463912	0.01175537509306583
+45	0.1888447419002717	0.7937432801187444	0.01741197798098377
+46	0.003343870943690294	0.5688485365439184	0.4278075925123913
+47	0.6298809126821926	0.351819105846624	0.01829998147118346
+48	0.02194748498208588	0.9528719644963803	0.02518055052153371
+49	0.0134809843973319	0.9158104804105686	0.07070853519209956
+50	0.2236040354818131	0.7658234826035886	0.01057248191459843
+51	2.468859571725732e-08	0.1149297558177613	0.8850702194936431
+52	0.8822943389422144	0.1174844188587244	0.0002212421990611336
+53	0.8212069491143118	0.1780292921423463	0.0007637587433421098
+54	0.01093910360053978	0.8913143071842399	0.09774658921522028
+55	0.282805088273554	0.7108876480767123	0.006307263649733691
+56	0.9056028772991068	0.09431812117689807	7.90015239951085e-05
+57	0.007004517863302363	0.7889401998535196	0.2040552822831782
+58	0.8240969291044774	0.1751797693467425	0.0007233015487799649
+59	0.004557214888219164	0.6654368937387614	0.3300058913730194
+60	0.5894432210638682	0.3812771655644528	0.02927961337167887
+61	0.1074604023879576	0.8109326188920644	0.08160697871997806
+62	0.06840301972803904	0.7387872167497551	0.1928097635222059
+63	0.01235430274922097	0.9108749443132855	0.07677075293749337
+64	0.02378845196344652	0.9548860769347108	0.02132547110184272
+65	0.02140041410781593	0.9513261243422556	0.02727346154992843
+66	2.652470045730108e-12	0.001536740791252874	0.9984632592060946
+67	0.4293436714667264	0.4449305495588506	0.125725778974423
+68	3.062451890924269e-05	0.3963884929404985	0.6035808825405923
+69	4.203558676524377e-07	0.1678273598150428	0.8321722198290895
+70	0.1285094505873794	0.8164374278781318	0.05505312153448871
+71	0.08884834523452112	0.788938136071633	0.1222135186938459
+72	0.02485971214175629	0.9541684094196913	0.02097187843855246
+73	0.0528392461971367	0.6705724187016827	0.2765883351011806
+74	0.8455383320112253	0.1540050744936643	0.0004565934951104923
+75	0.004956233951116152	0.6891677512204614	0.3058760148284225
+76	1.264448454146694e-06	0.3261147865980806	0.6738839489534654
+77	0.8884301422405665	0.1113406774536455	0.0002291803057879414
+78	0.6483585631910426	0.3367108876687838	0.01493054914017372
+79	0.3317383031925827	0.6650554670575105	0.00320622974990686
+80	0.8875671919014441	0.1122826168668424	0.0001501912317134283
+81	0.719199494586421	0.2756481086936792	0.005152396719899931
+82	0.6317584308515358	0.3505455769133323	0.01769599223513166
+83	9.381919643271678e-07	0.1870973766348126	0.8129016851732231
+84	9.256665300775749e-05	0.7230123512024462	0.2768950821445459
+85	0.8808332302148698	0.1189714722975506	0.0001952974875795364
+86	0.8822943389422144	0.1174844188587244	0.0002212421990611336
+87	1.493406683878595e-08	0.01448722859192672	0.9855127564740064
+88	0.6317584308515358	0.3505455769133323	0.01769599223513166
+89	0.811740904859526	0.1870917845442495	0.001167310596224406
+90	0.007439120532096984	0.8070742412353221	0.1854866382325809
+91	0.475629118295347	0.4364206143110693	0.08795026739358358
diff --git a/tests/hbltests/Results/HIVSweden.out_MODEL_-1.nex b/tests/hbltests/Results/HIVSweden.out_MODEL_-1.nex
index c6f922d03..345ebbd77 100644
--- a/tests/hbltests/Results/HIVSweden.out_MODEL_-1.nex
+++ b/tests/hbltests/Results/HIVSweden.out_MODEL_-1.nex
@@ -37,8 +37,8 @@ END;
 BEGIN HYPHY;
 
 
-global c=0.9053687405962891;
-global kappa=0.4055902889768502;
+global c=0.9053687083867886;
+global kappa=0.4055902861703429;
 modelMatrix={61,61};
 modelMatrix[0][1]:=kappa*c*t;
 modelMatrix[0][2]:=t;
@@ -638,29 +638,29 @@ ACCEPT_ROOTED_TREES=0;
 UseModel (theModel);
 Tree givenTree=(((U68496,U68497)Node3,(U68501,(U68502,(U68503,((((U68504,U68506)Node15,U68505)Node14,U68507)Node13,U68508)Node12)Node10)Node8)Node6)Node2,U68498,(U68499,U68500)Node6_1);
 
-givenTree.U68496.t=0.1932772918590503;
-givenTree.U68497.t=0.643982465056076;
-givenTree.Node3.t=1.242380162273858;
-givenTree.U68501.t=1.070573211757976;
-givenTree.U68502.t=1.555699524211797;
-givenTree.U68503.t=1.572551570577959;
-givenTree.U68504.t=0.1543465830348936;
-givenTree.U68506.t=0.5722610548396673;
-givenTree.Node15.t=0.4521342543764549;
-givenTree.U68505.t=0.2463547921390629;
-givenTree.Node14.t=0.2981011795265192;
-givenTree.U68507.t=0.4971516409620237;
-givenTree.Node13.t=0.2830238149773393;
-givenTree.U68508.t=1.612939402336149;
-givenTree.Node12.t=0.4111635119993902;
-givenTree.Node10.t=0.5204945671868649;
-givenTree.Node8.t=0.2739113550785894;
-givenTree.Node6.t=0.39109971254045;
-givenTree.Node2.t=0.192122168657457;
-givenTree.U68498.t=0.4097794060919068;
-givenTree.U68499.t=0.2309458844863399;
-givenTree.U68500.t=0.8122432466308773;
-givenTree.Node6_1.t=0.6787883745417747;
+givenTree.U68496.t=0.1932772963745281;
+givenTree.U68497.t=0.6439824826960953;
+givenTree.Node3.t=1.242380188093008;
+givenTree.U68501.t=1.070573236061633;
+givenTree.U68502.t=1.555699559158634;
+givenTree.U68503.t=1.572551629409277;
+givenTree.U68504.t=0.1543465863405398;
+givenTree.U68506.t=0.5722610718109247;
+givenTree.Node15.t=0.4521342669889602;
+givenTree.U68505.t=0.2463547989083778;
+givenTree.Node14.t=0.2981011966952727;
+givenTree.U68507.t=0.4971516512045175;
+givenTree.Node13.t=0.2830238933165307;
+givenTree.U68508.t=1.61293943451744;
+givenTree.Node12.t=0.4111635089345889;
+givenTree.Node10.t=0.5204945721494226;
+givenTree.Node8.t=0.2739113713935601;
+givenTree.Node6.t=0.39109971970209;
+givenTree.Node2.t=0.1921221738456284;
+givenTree.U68498.t=0.4097794151846853;
+givenTree.U68499.t=0.2309458913872281;
+givenTree.U68500.t=0.8122432676809466;
+givenTree.Node6_1.t=0.6787883949859681;
 DataSet ds = ReadDataFile(USE_NEXUS_FILE_DATA);
 DataSetFilter filteredData = CreateFilter(ds,3,"0-272","0,1,5-8,10,9,11,12,2-4","TAA,TAG,TGA");
 ASSUME_REVERSIBLE_MODELS=0;
diff --git a/tests/hbltests/Results/HIVSweden.out_MODEL_0.nex b/tests/hbltests/Results/HIVSweden.out_MODEL_0.nex
index 5c99a6f17..d4f9f60cb 100644
--- a/tests/hbltests/Results/HIVSweden.out_MODEL_0.nex
+++ b/tests/hbltests/Results/HIVSweden.out_MODEL_0.nex
@@ -37,11 +37,11 @@ END;
 BEGIN HYPHY;
 
 
-global W=0.07841413289607745;
+global W=0.07843652795974788;
 W:<1;
-global P=0.4833143321109446;
+global P=0.4829453153298657;
 P:<1;
-global kappa=0.3859671725353098;
+global kappa=0.3864165111352775;
 
 c.weights={1,2};
 c.weights[0][0]:=P;
@@ -654,29 +654,29 @@ ACCEPT_ROOTED_TREES=0;
 UseModel (theModel);
 Tree givenTree=(((U68496,U68497)Node3,(U68501,(U68502,(U68503,((((U68504,U68506)Node15,U68505)Node14,U68507)Node13,U68508)Node12)Node10)Node8)Node6)Node2,U68498,(U68499,U68500)Node6_1);
 
-givenTree.U68496.t=0.3019868795452008;
-givenTree.U68497.t=1.060303689299438;
-givenTree.Node3.t=2.087086139102;
-givenTree.U68501.t=1.76765256783072;
-givenTree.U68502.t=2.680492507749348;
-givenTree.U68503.t=2.702229081386573;
-givenTree.U68504.t=0.2432357442362262;
-givenTree.U68506.t=0.9044656237731131;
-givenTree.Node15.t=0.7158520457122713;
-givenTree.U68505.t=0.3957103842042503;
-givenTree.Node14.t=0.5004222842328082;
-givenTree.U68507.t=0.7752220423105534;
-givenTree.Node13.t=0.2759340027399538;
-givenTree.U68508.t=2.851296690424618;
-givenTree.Node12.t=0.7487808219465917;
-givenTree.Node10.t=0.8231033441506819;
-givenTree.Node8.t=0.3088428014417186;
-givenTree.Node6.t=0.6411477022454601;
-givenTree.Node2.t=0.3046465337039452;
-givenTree.U68498.t=0.6390986266843887;
-givenTree.U68499.t=0.3598159761956471;
-givenTree.U68500.t=1.320092572770488;
-givenTree.Node6_1.t=1.137305847275062;
+givenTree.U68496.t=0.3016777767498887;
+givenTree.U68497.t=1.058409360343901;
+givenTree.Node3.t=2.08013961029416;
+givenTree.U68501.t=1.764267438750952;
+givenTree.U68502.t=2.670061670837944;
+givenTree.U68503.t=2.700870351825198;
+givenTree.U68504.t=0.2430764173949833;
+givenTree.U68506.t=0.9034886180987204;
+givenTree.Node15.t=0.7157098400780856;
+givenTree.U68505.t=0.3954395994916897;
+givenTree.Node14.t=0.4993961571255977;
+givenTree.U68507.t=0.7736874142613375;
+givenTree.Node13.t=0.276686065842527;
+givenTree.U68508.t=2.844041601649456;
+givenTree.Node12.t=0.7474772558959789;
+givenTree.Node10.t=0.8233021855120572;
+givenTree.Node8.t=0.3100594201037387;
+givenTree.Node6.t=0.6408641029757065;
+givenTree.Node2.t=0.3046545921548557;
+givenTree.U68498.t=0.6387102496615624;
+givenTree.U68499.t=0.3595616425563687;
+givenTree.U68500.t=1.31867988076752;
+givenTree.Node6_1.t=1.137365801976078;
 DataSet ds = ReadDataFile(USE_NEXUS_FILE_DATA);
 DataSetFilter filteredData = CreateFilter(ds,3,"0-272","0,1,5-8,10,9,11,12,2-4","TAA,TAG,TGA");
 ASSUME_REVERSIBLE_MODELS=0;
diff --git a/tests/hbltests/Results/HIVSweden.out_MODEL_1.nex b/tests/hbltests/Results/HIVSweden.out_MODEL_1.nex
index 0fdccdf78..eb448e519 100644
--- a/tests/hbltests/Results/HIVSweden.out_MODEL_1.nex
+++ b/tests/hbltests/Results/HIVSweden.out_MODEL_1.nex
@@ -37,15 +37,15 @@ END;
 BEGIN HYPHY;
 
 
-global P2=0.7086571320382176;
+global P2=0.7091890882964867;
 P2:<1;
-global W_1=0.05881756131775031;
+global W_1=0.05836173276024324;
 W_1:<1;
-global W_2=3.635394261991776;
+global W_2=3.637969705855394;
 W_2:>1;
-global P1=0.3752208738212064;
+global P1=0.3746328183428458;
 P1:<1;
-global kappa=0.3595859416144088;
+global kappa=0.359547071282534;
 
 c.weights={1,3};
 c.weights[0][0]:=P1;
@@ -660,29 +660,29 @@ ACCEPT_ROOTED_TREES=0;
 UseModel (theModel);
 Tree givenTree=(((U68496,U68497)Node3,(U68501,(U68502,(U68503,((((U68504,U68506)Node15,U68505)Node14,U68507)Node13,U68508)Node12)Node10)Node8)Node6)Node2,U68498,(U68499,U68500)Node6_1);
 
-givenTree.U68496.t=0.1821192790979838;
-givenTree.U68497.t=0.6321394500123422;
-givenTree.Node3.t=1.245080511201583;
-givenTree.U68501.t=1.054344150733757;
-givenTree.U68502.t=1.750026688242964;
-givenTree.U68503.t=1.776229064972074;
-givenTree.U68504.t=0.1466190694742395;
-givenTree.U68506.t=0.5212616525946417;
-givenTree.Node15.t=0.4176673911941267;
-givenTree.U68505.t=0.241487809690917;
-givenTree.Node14.t=0.2760969899208042;
-givenTree.U68507.t=0.4779878429167616;
-givenTree.Node13.t=0.005926806848416143;
-givenTree.U68508.t=1.943186033352408;
-givenTree.Node12.t=0.5089108795301405;
-givenTree.Node10.t=0.5085100270232016;
-givenTree.Node8.t=0.1808510348041232;
-givenTree.Node6.t=0.3945987177203757;
-givenTree.Node2.t=0.1899356130052154;
-givenTree.U68498.t=0.3343574439167225;
-givenTree.U68499.t=0.1881468418422207;
-givenTree.U68500.t=0.8010667515952579;
-givenTree.Node6_1.t=0.7504188627354821;
+givenTree.U68496.t=0.1821113220484527;
+givenTree.U68497.t=0.6316817951793845;
+givenTree.Node3.t=1.246078874282207;
+givenTree.U68501.t=1.054399281533476;
+givenTree.U68502.t=1.749683212917501;
+givenTree.U68503.t=1.776615432547105;
+givenTree.U68504.t=0.1466531311564719;
+givenTree.U68506.t=0.5208909511938736;
+givenTree.Node15.t=0.4175019163743378;
+givenTree.U68505.t=0.2416415468045435;
+givenTree.Node14.t=0.2760241601189965;
+givenTree.U68507.t=0.4780632082426846;
+givenTree.Node13.t=0.006038075586898496;
+givenTree.U68508.t=1.943589922192904;
+givenTree.Node12.t=0.5088019660443692;
+givenTree.Node10.t=0.5079730105635974;
+givenTree.Node8.t=0.1809257643034209;
+givenTree.Node6.t=0.3946675503522505;
+givenTree.Node2.t=0.1897537804085054;
+givenTree.U68498.t=0.3347160083706835;
+givenTree.U68499.t=0.1882826930766557;
+givenTree.U68500.t=0.8008037461971568;
+givenTree.Node6_1.t=0.7498993417480263;
 DataSet ds = ReadDataFile(USE_NEXUS_FILE_DATA);
 DataSetFilter filteredData = CreateFilter(ds,3,"0-272","0,1,5-8,10,9,11,12,2-4","TAA,TAG,TGA");
 ASSUME_REVERSIBLE_MODELS=0;
diff --git a/tests/hbltests/UnitTests/HBLCommands/DataSetFilter.bf b/tests/hbltests/UnitTests/HBLCommands/DataSetFilter.bf
index 67fea1570..b8a9eebe5 100644
--- a/tests/hbltests/UnitTests/HBLCommands/DataSetFilter.bf
+++ b/tests/hbltests/UnitTests/HBLCommands/DataSetFilter.bf
@@ -4,13 +4,13 @@ runATest ();
 
 function getTestName () {
   return "DataSetFilter";
-}		
+}
 
 
 function runTest () {
 	ASSERTION_BEHAVIOR = 1; /* print warning to console and go to the end of the execution list */
-	testResult = TRUE;
-  
+	testResult = 0;
+
 
   //---------------------------------------------------------------------------------------------------------
   // SIMPLE FUNCTIONALITY
@@ -20,7 +20,7 @@ function runTest () {
 
   DataSet cd2nex = ReadDataFile (PATH_TO_CURRENT_BF + '/../../data/CD2.nex');
 
-  DataSetFilter unChanged = CreateFilter (cd2nex,1,"","","");  
+  DataSetFilter unChanged = CreateFilter (cd2nex,1,"","","");
   DataSetFilter removedAAA = CreateFilter (cd2nex,3,"","","AAA");
   DataSetFilter onlyFiveThroughTen = CreateFilter (cd2nex,1,"4-9");
   DataSetFilter onlyLiveStock = CreateFilter (cd2nex,1,"","4-6");
@@ -31,20 +31,20 @@ function runTest () {
   HarvestFrequencies (freqsOnlyFiveThroughTen, onlyFiveThroughTen, 1, 1, 1);
   HarvestFrequencies (freqsOnlyLiveStock, onlyLiveStock, 1, 1, 1);
 
-  assert(Abs((freqsUnChanged[0] - freqsRemovedAAA[0]) - 0.0392668) < 0.0001 , "Failed to remove AAA codons with DataSetFilter");
+  assert(Abs((freqsUnChanged[0] - freqsRemovedAAA[0]) - 0.073822) < 0.0001 , "Failed to remove AAA codons with DataSetFilter");
   assert(Abs((freqsUnChanged[0] - freqsOnlyFiveThroughTen[0]) + 0.0145053) < 0.001, "Failed to filter based on site index with DataSetFilter");
   assert(Abs((freqsUnChanged[0] - freqsOnlyLiveStock[0]) - 0.00130718) < 0.001 , "Failed to filter based on sequence index with DataSetFilter");
-  
+
   //---------------------------------------------------------------------------------------------------------
   // ERROR HANDLING
   //---------------------------------------------------------------------------------------------------------
   list1 = {'key1':'val1', 'key2':'val2'};
-  
+
   assert (runCommandWithSoftErrors ("DataSetFilter test1 = CreateFilter (list1,1,'','4-6');", "Expected a non-empty string matrix as the FILTER_NAMES argument in call to CreateFilter"), "Failed error checking for trying to perform a dataset filter on a non-file object");
-  
+
   //TODO{low-priority}: The line below results in "Floating point exception: 8"
   //assert (runCommandWithSoftErrors ("DataSetFilter test2 = CreateFilter (cd2nex,list1,'','4-6');", "was expected to be a numerical argument"), "Failed error checking for trying to perform a dataset filter on a non-file object");
-  
+
   testResult = 1;
 
   return testResult;
diff --git a/tests/hbltests/UnitTests/HBLCommands/DeleteObject.bf b/tests/hbltests/UnitTests/HBLCommands/DeleteObject.bf
index 30148d0f0..75f7cd774 100644
--- a/tests/hbltests/UnitTests/HBLCommands/DeleteObject.bf
+++ b/tests/hbltests/UnitTests/HBLCommands/DeleteObject.bf
@@ -15,16 +15,14 @@ function getTestedFunctions ()
 
 function runTest () {
 	ASSERTION_BEHAVIOR = 1; /* print warning to console and go to the end of the execution list */
-	testResult  	   = 1;
+	testResult  	   = 0;
 
-    ExecuteAFile (PATH_TO_CURRENT_BF  + "res" + DIRECTORY_SEPARATOR + "test_likefunc2.nex");
-    ExecuteAFile (PATH_TO_CURRENT_BF  + "res" + DIRECTORY_SEPARATOR + "test_likefunc.nex", {}, "a_prefix");
-    ExecuteAFile (PATH_TO_CURRENT_BF  + "res" + DIRECTORY_SEPARATOR + "test_likefunc.nex");
+  ExecuteAFile (PATH_TO_CURRENT_BF  + "res" + DIRECTORY_SEPARATOR + "test_likefunc2.nex");
+  ExecuteAFile (PATH_TO_CURRENT_BF  + "res" + DIRECTORY_SEPARATOR + "test_likefunc.nex", {}, "a_prefix");
+  ExecuteAFile (PATH_TO_CURRENT_BF  + "res" + DIRECTORY_SEPARATOR + "test_likefunc.nex");
 	testResult  	   = TRUE;
 
 
 
 	return 1;
 }
-
-
diff --git a/tests/hbltests/UnitTests/HBLCommands/GetDataInfo.bf b/tests/hbltests/UnitTests/HBLCommands/GetDataInfo.bf
index d024b4efe..2b5d1f430 100644
--- a/tests/hbltests/UnitTests/HBLCommands/GetDataInfo.bf
+++ b/tests/hbltests/UnitTests/HBLCommands/GetDataInfo.bf
@@ -1,10 +1,10 @@
 function getTestName () {
 	return "GetDataInfo";
-}		
+}
 
 function getTestedFunctions () {
 	return {{"_ExecutionList::ExecuteCase46","_DataSetFilter::FindUniqueSequences"}};
-}	
+}
 
 function runTest () {
 	/* define some auxiliary objects here */
@@ -14,9 +14,9 @@ function runTest () {
 
 	DataSet 			fiveSeqs = ReadDataFile (PATH_TO_CURRENT_BF + "../../data/5.fas");	
 	DataSetFilter		nucF	 = CreateFilter (fiveSeqs,1);
-	
+
 	GetDataInfo 		(seqInfo, nucF, -1);
-		
+
 	assert (seqInfo["UNIQUE_SEQUENCES"] == 5, "Expected 5 unique sequences with strict filtering");
 	GetDataInfo 		(seqInfo, nucF, -2);
 	assert (seqInfo["UNIQUE_SEQUENCES"] == 4, "Expected 4 unique sequences with strict+gap filtering");
@@ -26,11 +26,11 @@ function runTest () {
 	assert (seqInfo["UNIQUE_SEQUENCES"] == 2, "Expected 2 unique sequences with partial match filtering");
 
 	DataSetFilter		dinucF	 = CreateFilter (fiveSeqs,2);
-		
+
 	GetDataInfo 		(seqInfo, dinucF, -2);
 	assert (seqInfo["UNIQUE_SEQUENCES"] == 5, "Expected 5 unique sequences with strict+gap filtering (dinuc)");
-		
-	testResult = 1;	
+
+	testResult = 1;
 	return testResult;
 }
 
diff --git a/tests/hbltests/UnitTests/HBLCommands/ReplicateConstraint.bf b/tests/hbltests/UnitTests/HBLCommands/ReplicateConstraint.bf
index ae1541dc2..53da98319 100644
--- a/tests/hbltests/UnitTests/HBLCommands/ReplicateConstraint.bf
+++ b/tests/hbltests/UnitTests/HBLCommands/ReplicateConstraint.bf
@@ -171,7 +171,6 @@ function runTest () {
 		   );
 
 
-
 	testResult = 1;
 
 	return testResult;
diff --git a/tests/hbltests/libv3/ABSREL.wbf b/tests/hbltests/libv3/ABSREL.wbf
index eb70f97e2..f1131f76d 100644
--- a/tests/hbltests/libv3/ABSREL.wbf
+++ b/tests/hbltests/libv3/ABSREL.wbf
@@ -8,13 +8,6 @@ if (+version >= 2.4) {
 }
 LoadFunctionLibrary ("shared.bf");
 
-/*
-LoadFunctionLibrary ("libv3/IOFunctions.bf");
-fscanf ("data/CD2.nex.ABSREL.json","Raw",json);
-absrel.json = Eval (json);
-*/
-
-
 assert (check_value (
     ((absrel.json["fits"])["Full adaptive model"])["Log Likelihood"], -3415.02, 0.001), "Incorrect log-likelihood for the full adaptive model");
 
@@ -23,13 +16,14 @@ assert (check_value (
 
 assert (check_value (
     ((absrel.json["test results"])["tested"]),6, 0.001), "Incorrect number of total tests");
-    
 
 
-test.expected_positives = utility.MatrixToDict({{"Node2","Node8","BEAVIS"}});
+
+test.expected_positives = utility.MatrixToDict({{"Node2","Node8"}});
 test.lrts = 0;
 
 function confirm_branch (branch, p, dict) {
+
     if (p == None) {
         p = 1;
     }
@@ -47,15 +41,12 @@ function confirm_branch (branch, p, dict) {
     return false;
 }
 
-utility.ForEachPair ((absrel.json["branch attributes"])["0"],"_key_", "_value_",
-                 "  
-                    if (confirm_branch (_key_, _value_['Corrected P-value'], test.expected_positives)) {
-                    test.lrts  += _value_['LRT'];
-                   }
-                ");
-
+for (_key_, _value_; in; (absrel.json["branch attributes"])["0"]) {
+  if (_value_/'Corrected P-value') {
+    if (confirm_branch (_key_, _value_['Corrected P-value'], test.expected_positives)) {
+      test.lrts  += _value_['LRT'];
+    }
+  }
+}
 assert (check_value (
     test.lrts, 15.67, 0.05), "More than 5% difference in cumulative LRT for positively selected sites");
-    
-
-