diff --git a/CMakeLists.txt b/CMakeLists.txt
index b170af57d..d41a1b7b4 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -579,7 +579,7 @@ add_test (NAME SLAC COMMAND HYPHYMP tests/hbltests/libv3/SLAC.wbf)
 add_test (NAME SLAC-PARTITIONED COMMAND HYPHYMP tests/hbltests/libv3/SLAC-partitioned.wbf)
 add_test (NAME FEL COMMAND HYPHYMP tests/hbltests/libv3/FEL.wbf)
 add_test (MEME HYPHYMP tests/hbltests/libv3/MEME.wbf)
-add_test (MEME-PARTITIONED HYPHYMP tests/hbltests/libv3/MEME.wbf)
+add_test (MEME-PARTITIONED HYPHYMP tests/hbltests/libv3/MEME-partitioned.wbf)
 add_test (BUSTED HYPHYMP tests/hbltests/libv3/BUSTED.wbf)
 add_test (BUSTED-SRV HYPHYMP tests/hbltests/libv3/BUSTED-SRV.wbf)
 add_test (RELAX HYPHYMP tests/hbltests/libv3/RELAX.wbf)
diff --git a/res/TemplateBatchFiles/SelectionAnalyses/FEL.bf b/res/TemplateBatchFiles/SelectionAnalyses/FEL.bf
index a008bd8cb..bf13802ce 100644
--- a/res/TemplateBatchFiles/SelectionAnalyses/FEL.bf
+++ b/res/TemplateBatchFiles/SelectionAnalyses/FEL.bf
@@ -339,7 +339,7 @@ lfunction fel.handle_a_site (lf, filter_data, partition_index, pattern_info, mod
     ^"fel.alpha_scaler" = (^"fel.alpha_scaler" + 3*^"fel.beta_scaler_test")/4;
     parameters.SetConstraint ("fel.beta_scaler_test","fel.alpha_scaler", "");
 
-    //Optimize (results, ^lf, {"OPTIMIZATION_METHOD" : "nedler-mead", OPTIMIZATION_PRECISION: 1e-5});
+    //Optimize (results, ^lf, {"OPTIMIZATION_METHOD" : "nedler-mead"});
     Optimize (results, ^lf);
 
     Null = estimators.ExtractMLEs (lf, model_mapping);
diff --git a/res/TemplateBatchFiles/SelectionAnalyses/MEME.bf b/res/TemplateBatchFiles/SelectionAnalyses/MEME.bf
index e6df92340..a56fd8e42 100644
--- a/res/TemplateBatchFiles/SelectionAnalyses/MEME.bf
+++ b/res/TemplateBatchFiles/SelectionAnalyses/MEME.bf
@@ -118,7 +118,7 @@ This table is meant for HTML rendering in the results web-app; can use HTML char
 is 'pop-over' explanation of terms. This is ONLY saved to the JSON file. For Markdown screen output see
 the next set of variables.
 */
-meme.table_screen_output  = {{"Codon", "Partition", "alpha", "beta+", "p+", "LRT", "Episodic selection detected?", "# branches"}};
+meme.table_screen_output  = {{"Codon", "Partition", "alpha", "beta+", "p+", "LRT", "Episodic selection detected?", "# branches", "Most common codon substitutions at this site"}};
 meme.table_output_options = {terms.table_options.header : TRUE, terms.table_options.minimum_column_width: 12, terms.table_options.align : "center"};
 
 
@@ -279,6 +279,8 @@ parameters.SetRange ("meme.site_mixture_weight", terms.range_almost_01);
 
 meme.report.count = {{0}};
 
+meme.site.composition.string = "";
+
 meme.report.positive_site = {{"" + (1+((meme.filter_specification[meme.report.partition])[terms.data.coverage])[meme.report.site]),
                                     meme.report.partition + 1,
                                     Format(meme.report.row[0],7,3),
@@ -286,7 +288,8 @@ meme.report.positive_site = {{"" + (1+((meme.filter_specification[meme.report.pa
                                     Format(meme.report.row[4],7,3),
                                     Format(meme.report.row[5],7,3),
                                     "Yes, p = " + Format(meme.report.row[6],7,4),
-                                    Format(meme.report.row[7],0,0)
+                                    Format(meme.report.row[7],0,0),
+                                    meme.site.composition.string
 }};
 
 meme.site_results = {};
@@ -575,13 +578,8 @@ lfunction meme.handle_a_site (lf_fel, lf_bsrel, filter_data, partition_index, pa
     ancestral_info = ancestral.build (lf_bsrel,0,FALSE);
 
     //TODO
-    branch_substitution_information = selection.substitution_mapper (ancestral_info ["MATRIX"],
-                                                      ancestral_info ["TREE_AVL"],
-                                                      ancestral_info ["AMBIGS"],
-                                                      ^"meme.pairwise_counts",
-                                                      ancestral_info ["MAPPING"],
-                                                      (^"meme.codon_data_info")[utility.getGlobalValue("terms.code")]);
-
+    branch_substitution_information = (ancestral.ComputeSubstitutionBySite (ancestral_info,0,None))[^"terms.substitutions"];
+    
 
     DeleteObject (ancestral_info);
 
@@ -607,7 +605,7 @@ lfunction meme.handle_a_site (lf_fel, lf_bsrel, filter_data, partition_index, pa
         LFCompute (^lf_bsrel,LF_DONE_COMPUTE);
 
         ^"meme.site_beta_plus" := ^"meme.site_alpha";
-        Optimize (results, ^lf_bsrel);
+        Optimize (results, ^lf_bsrel, {"OPTIMIZATION_METHOD" : "nedler-mead"});
         //io.SpoolLF (lf_bsrel, "/tmp/meme.debug", "MEME-null");
 
         Null = estimators.ExtractMLEs (lf_bsrel, model_mapping);
@@ -666,6 +664,60 @@ function meme.report.echo (meme.report.site, meme.report.partition, meme.report.
 
 lfunction meme.store_results (node, result, arguments) {
 
+    sub_profile = result[utility.getGlobalValue("terms.branch_selection_attributes")];
+    
+/**
+
+{
+ "AGA":{
+   "AGG":{
+     "0":"Node1"
+    }
+  },
+ "AGG":{
+   "AAA":{
+     "0":"Node3",
+     "1":"HORSE"
+    },
+   "ACG":{
+     "0":"Node12"
+    },
+   "GCA":{
+     "0":"CAT"
+    }
+  }
+}
+
+**/
+
+    if (None != sub_profile) {
+        total_sub_count = 0;
+    
+        sub_by_count = {};
+        for (i, v; in; sub_profile) {
+            for (j, b; in; v) {
+                c = Abs (b);
+                total_sub_count += c;
+                if ((sub_by_count/c)==0) {
+                    sub_by_count [c] = {};
+                }
+                sub_by_count[c] + (i + ">" + j);
+            }
+        }
+        
+        sorted_subs = {Abs (sub_by_count), 1};
+        j = 0;
+        for (i, v; in; sub_by_count) {
+            sorted_subs[j] = -(+i);
+            j += 1; 
+        }
+        sub_profile = {};
+        for (i; in; sorted_subs % 0) {
+            sub_profile + ("[" + (-i) + "]" + Join (",",sub_by_count [-i]));
+        }
+        ^"meme.site.composition.string" = Join ("|", sub_profile);
+    }
+
     partition_index = arguments [3];
     pattern_info    = arguments [4];
 
@@ -730,6 +782,7 @@ lfunction meme.store_results (node, result, arguments) {
     sites_mapping_to_pattern = pattern_info[utility.getGlobalValue("terms.data.sites")];
     sites_mapping_to_pattern.count = utility.Array1D (sites_mapping_to_pattern);
 
+
     for (i = 0; i < sites_mapping_to_pattern.count; i+=1) {
         site_index = sites_mapping_to_pattern[i];
         ((^"meme.site_results")[partition_index])[site_index] = result_row;
@@ -752,3 +805,4 @@ lfunction meme.store_results (node, result, arguments) {
 
     //assert (0);
 }
+
diff --git a/res/TemplateBatchFiles/SelectionAnalyses/SLAC.bf b/res/TemplateBatchFiles/SelectionAnalyses/SLAC.bf
index 2ed782984..89e5425b3 100644
--- a/res/TemplateBatchFiles/SelectionAnalyses/SLAC.bf
+++ b/res/TemplateBatchFiles/SelectionAnalyses/SLAC.bf
@@ -246,6 +246,7 @@ for (slac.i = 0; slac.i < Abs (slac.filter_specification); slac.i += 1) {
     slac.table_output_options[terms.table_options.header] = TRUE;
 
     slac.ancestors         = ancestral.build (slac.partitioned_mg_results[terms.likelihood_function], slac.i, None);
+    
     slac.results           [slac.i] = slac.compute_the_counts (slac.ancestors["MATRIX"], slac.ancestors["TREE_AVL"], slac.ancestors["AMBIGS"], slac.selected_branches[slac.i], slac.counts);
 
     slac.partition_sites   = utility.Array1D ((slac.filter_specification[slac.i])[terms.data.coverage]);
@@ -575,11 +576,11 @@ lfunction slac.compute_the_counts (matrix, tree, lookup, selected_branches, coun
             relative_branch_length = selected_branches_lengths[i] / selected_branch_total_length;
 
             parent_branch          = selected_branches_parents[i];
-
+ 
              for (s = 0; s < site_count; s += 1) {
                 this_state      = matrix[this_branch][s];
                 parent_state    = matrix[parent_branch][s];
-
+                
                 // parent state can only be resolved or --- (-1)
 
                 if (this_state >= 0) { // child fully resolved (this means that the parent is fully resolved as well)
diff --git a/res/TemplateBatchFiles/SelectionAnalyses/contrast-fel.bf b/res/TemplateBatchFiles/SelectionAnalyses/contrast-fel.bf
index 0e7737472..00f09c5f7 100644
--- a/res/TemplateBatchFiles/SelectionAnalyses/contrast-fel.bf
+++ b/res/TemplateBatchFiles/SelectionAnalyses/contrast-fel.bf
@@ -748,7 +748,7 @@ lfunction fel.handle_a_site (lf, filter_data, partition_index, pattern_info, mod
         );
     }
 
-    Optimize (results, ^lf);
+    Optimize (results, ^lf, {"OPTIMIZATION_METHOD" : "nedler-mead"});
     Null = estimators.ExtractMLEs (lf, model_mapping);
     Null [utility.getGlobalValue("terms.fit.log_likelihood")] = results[1][0];
 
@@ -761,7 +761,7 @@ lfunction fel.handle_a_site (lf, filter_data, partition_index, pattern_info, mod
 
                 estimators.RestoreLFStateFromSnapshot (lf_id, snapshot);
                 parameters.SetConstraint ((^"fel.scaler_parameter_names")[v1n],(^"fel.scaler_parameter_names")[v2n], "");
-                Optimize (results, ^lf);
+                Optimize (results, ^lf, {"OPTIMIZATION_METHOD" : "nedler-mead"});
                 pairwise[v1n + "|" + v2n] = estimators.ExtractMLEs (lf, model_mapping);
                 (pairwise[v1n + "|" + v2n])[utility.getGlobalValue("terms.fit.log_likelihood")] = results[1][0];
            }
diff --git a/res/TemplateBatchFiles/SelectionAnalyses/modules/selection_lib.ibf b/res/TemplateBatchFiles/SelectionAnalyses/modules/selection_lib.ibf
index cf84f98e5..1dddc9ceb 100644
--- a/res/TemplateBatchFiles/SelectionAnalyses/modules/selection_lib.ibf
+++ b/res/TemplateBatchFiles/SelectionAnalyses/modules/selection_lib.ibf
@@ -97,7 +97,7 @@ lfunction selection.substitution_mapper (matrix, tree, ambig_lookup, pairwise_co
         
         
         for (key, value; in; branch_info) {
-            result[key][bname] = value;
+            (result[key])[bname] = value;
         }
 
         //utility.ForEach (utility.Keys (branch_info), "slac.substituton_mapper.key",
diff --git a/res/TemplateBatchFiles/libv3/UtilityFunctions.bf b/res/TemplateBatchFiles/libv3/UtilityFunctions.bf
index a1a63591c..b46838392 100644
--- a/res/TemplateBatchFiles/libv3/UtilityFunctions.bf
+++ b/res/TemplateBatchFiles/libv3/UtilityFunctions.bf
@@ -1054,26 +1054,27 @@ function utility.sortStrings (_theList) {
 	return _gb_sortedStrings;
 }
 
-/*
+
 function utility.FinishAndPrintProfile () {
-#profile _hyphy_profile_dump;
+    ExecuteCommands ('
+        #profile _hyphy_profile_dump;
 
-    stats  			= _hyphy_profile_dump["STATS"];
-    _profile_summer = ({1,Rows(stats)}["1"]) * stats;
-    _instructions   = _hyphy_profile_dump["INSTRUCTION"];
-    _indices	    = _hyphy_profile_dump["INSTRUCTION INDEX"];
+            stats  			= _hyphy_profile_dump["STATS"];
+            _profile_summer = ({1,Rows(stats)}["1"]) * stats;
+            _instructions   = _hyphy_profile_dump["INSTRUCTION"];
+            _indices	    = _hyphy_profile_dump["INSTRUCTION INDEX"];
 
-    fprintf (stdout, "\nTotal run time (seconds)      : ", Format(_profile_summer[1],15,6),
-                     "\nTotal number of steps         : ", Format(_profile_summer[0],15,0), "\n\n");
+            fprintf (stdout, "\nTotal run time (seconds)      : ", Format(_profile_summer[1],15,6),
+                             "\nTotal number of steps         : ", Format(_profile_summer[0],15,0), "\n\n");
 
-    to_sort        =  stats["-_MATRIX_ELEMENT_VALUE_*_MATRIX_ELEMENT_COLUMN_+(_MATRIX_ELEMENT_COLUMN_==0)*_MATRIX_ELEMENT_ROW_"] % 1;
+            to_sort        =  stats["-_MATRIX_ELEMENT_VALUE_*_MATRIX_ELEMENT_COLUMN_+(_MATRIX_ELEMENT_COLUMN_==0)*_MATRIX_ELEMENT_ROW_"] % 1;
 
-    for (k=0; k<Columns(_instructions); k=k+1)
-    {
-        k2 = to_sort[k][0];
-        fprintf (stdout, Format (_indices[k2],6,0), " : ", _instructions[k2], "\n\tCall count: ", stats[k2][0],
-                                                       "\n\tTime (seconds): ", stats[k2][1], "\n");
-    }
-    return 0;
+            for (k=0; k<Columns(_instructions); k=k+1)
+            {
+                k2 = to_sort[k][0];
+                fprintf (stdout, Format (_indices[k2],6,0), " : ", _instructions[k2], "\n\tCall count: ", stats[k2][0],
+                                                               "\n\tTime (seconds): ", stats[k2][1], "\n");
+            }
+    ');
 }
-*/
+
diff --git a/res/TemplateBatchFiles/libv3/all-terms.bf b/res/TemplateBatchFiles/libv3/all-terms.bf
index 4401f0b4a..3f7467db3 100644
--- a/res/TemplateBatchFiles/libv3/all-terms.bf
+++ b/res/TemplateBatchFiles/libv3/all-terms.bf
@@ -179,6 +179,7 @@ namespace terms{
         coverage       = "coverage";
         filter_string  = "filter-string";
         is_constant    = "is_constant";
+        characters     = "characters";
         pattern_id     = "pattern id";
         filename_to_index = "filename-to-index";
     }
diff --git a/res/TemplateBatchFiles/libv3/tasks/alignments.bf b/res/TemplateBatchFiles/libv3/tasks/alignments.bf
index 6dc9bd088..d9f02f6ff 100644
--- a/res/TemplateBatchFiles/libv3/tasks/alignments.bf
+++ b/res/TemplateBatchFiles/libv3/tasks/alignments.bf
@@ -453,6 +453,8 @@ lfunction alignments.CompressDuplicateSequences (filter_in, filter_out, rename)
     //DataSetFilter ^filter_out = CreateFilter (filter_in);
 }
 
+
+
 /**
  * Defines filters for multiple partitions
  * @name alignments.DefineFiltersForPartitions
@@ -796,6 +798,8 @@ lfunction alignments.Extract_site_patterns (data_filter) {
     site_characters = {};
     sequence_count = ^(data_filter + ".species");
     
+    GetDataInfo (filter_characters, ^data_filter, "CHARACTERS");
+    
     
     for (_site_index_, _pattern_; in; pattern_list) {
         utility.EnsureKey (site_info, _pattern_);
@@ -810,6 +814,10 @@ lfunction alignments.Extract_site_patterns (data_filter) {
                 }
             }
             site_characters = sc;
+            (site_info[_pattern_])[^"terms.data.characters"] = {};
+            for (i,v; in; site_characters) {
+                ((site_info[_pattern_])[^"terms.data.characters"])[filter_characters[+i]] = v;
+            }
             (site_info[_pattern_])[^"terms.data.is_constant"] = Abs (site_characters) <= 1;
         }
         
diff --git a/res/TemplateBatchFiles/libv3/tasks/ancestral.bf b/res/TemplateBatchFiles/libv3/tasks/ancestral.bf
index 950ad6b33..05069ec4f 100644
--- a/res/TemplateBatchFiles/libv3/tasks/ancestral.bf
+++ b/res/TemplateBatchFiles/libv3/tasks/ancestral.bf
@@ -267,7 +267,7 @@ lfunction ancestral._buildAncestralCacheInternal(_lfID, _lfComponentID, doSample
                         GetDataInfo(_bacCharState, _bacAF, _bacRowIndex - _bacFilterSequenceCount, _bacAncestralPatternMap[_bacSiteCounter]);
                     }
                     _bacResolutionCount = +_bacCharState;
-                    if (_bacResolutionCount == Columns(_bacCharHandles)) {
+                    if (_bacResolutionCount == Columns(_bacCharHandles) || _bacResolutionCount == 0) {
                             /* gap/full ambig */
                         if ( (reverse_mapping/(-1)) == FALSE) {
                             _bacHandledResolutions[_bacCurrentState] = -1;
@@ -538,10 +538,11 @@ lfunction ancestral.ComputeSubstitutionBySite (ancestral_data, site, branch_filt
         parent = (selected_branches[b])[1];
         own_state    = (ancestral_data["MATRIX"])[self][site];
         parent_state = (ancestral_data["MATRIX"])[parent][site];
-
+        
+        
         if  ((own_state != parent_state) && (own_state != -1) && (parent_state != -1)) {
-            own_state = (ancestral_data["CHARS"])[own_state];
-            parent_state = (ancestral_data["CHARS"])[parent_state];
+            own_state = (ancestral_data["MAPPING"])[own_state];
+            parent_state = (ancestral_data["MAPPING"])[parent_state];
             utility.EnsureKey (result, parent_state);
             utility.EnsureKey (result[parent_state], own_state);
             ((result[parent_state])[own_state]) + selected_branch_names[b];
diff --git a/src/core/alignment.cpp b/src/core/alignment.cpp
index 1ac9923a9..17c2358b1 100644
--- a/src/core/alignment.cpp
+++ b/src/core/alignment.cpp
@@ -994,6 +994,7 @@ double AlignStrings( char const * r_str
 
                     // if anything drops below 0, something bad happened
                     if ( i < 0 || j < 0 ) {
+                        delete [] edit_ops;
                         return -INFINITY;
                     }
 
diff --git a/src/core/batchlanruntime.cpp b/src/core/batchlanruntime.cpp
index e76b48edb..12d501493 100644
--- a/src/core/batchlanruntime.cpp
+++ b/src/core/batchlanruntime.cpp
@@ -541,13 +541,16 @@ bool      _ElementaryCommand::HandleGetDataInfo (_ExecutionList& current_program
 
                     if (site >=0 && site<filter_source->GetPatternCount()) {
                         if ( seq>=0 && seq<filter_source->NumberSpecies()) {
+                            bool count_gaps = hy_env::EnvVariableTrue(hy_env::harvest_frequencies_gap_options);
                             _Matrix             * res = new _Matrix (filter_source->GetDimension (true), 1, false, true);
 
                             bool                only_the_index = hy_env::EnvVariableTrue(hy_env::get_data_info_returns_only_the_index);
 
                             _String             character (filter_source->RetrieveState(site, seq));
-                            long                theValue = filter_source->Translate2Frequencies (character, res->theData,  true);
-
+                            
+                            
+                            long                theValue = filter_source->Translate2Frequencies (character, res->theData,  count_gaps);
+  
                             if (only_the_index) {
                                 receptacle->SetValue (new _Constant (theValue),false,true, NULL);
                                 DeleteObject     (res);
diff --git a/src/core/formula.cpp b/src/core/formula.cpp
index 7e5bd14dd..80cbd5947 100644
--- a/src/core/formula.cpp
+++ b/src/core/formula.cpp
@@ -2057,7 +2057,7 @@ bool _Formula::CheckSimpleTerm (HBLObjectRef thisObj) {
         if (oc != NUMBER) {
             if (oc ==MATRIX) {
                 _Matrix * mv = (_Matrix*)thisObj->Compute();
-                if (!mv->SparseDataStructure() && mv->IsIndependent ()) {
+                if (mv->is_dense() && mv->IsIndependent ()) {
                     return true;
                 }
             }
diff --git a/src/core/global_things.cpp b/src/core/global_things.cpp
index 9069ab1d4..dfdf062be 100644
--- a/src/core/global_things.cpp
+++ b/src/core/global_things.cpp
@@ -120,7 +120,7 @@ namespace hy_global {
                      kErrorStringDatasetRefIndexError ("Dataset index reference out of range"),
                      kErrorStringMatrixExportError    ("Export matrix called with a non-polynomial matrix argument"),
                      kErrorStringNullOperand          ("Attempting to operate on an undefined value; this is probably the result of an earlier 'soft' error condition"),
-                     kHyPhyVersion  = _String ("2.5.19"),
+                     kHyPhyVersion  = _String ("2.5.20"),
     
                     kNoneToken = "None",
                     kNullToken = "null",
diff --git a/src/core/include/matrix.h b/src/core/include/matrix.h
index d969f0e91..96f6e6fc6 100644
--- a/src/core/include/matrix.h
+++ b/src/core/include/matrix.h
@@ -142,7 +142,7 @@ class       _Matrix: public _MathObject {
     
     // matrix exp truncation precision
     
- 
+    bool        _validateCompressedStorage (void) const;
 
 public:
     
@@ -518,9 +518,6 @@ class       _Matrix: public _MathObject {
     long              MatrixType (void) {
         return storageType;
     }
-    bool              SparseDataStructure (void) {
-        return theIndex;
-    }
     
     template <typename CALLBACK, typename EXTRACTOR>  void ForEach (CALLBACK&& cbv, EXTRACTOR&& accessor) const {
         if (theIndex) {
diff --git a/src/core/likefunc.cpp b/src/core/likefunc.cpp
index 98309afa7..f846c0d52 100644
--- a/src/core/likefunc.cpp
+++ b/src/core/likefunc.cpp
@@ -5288,6 +5288,10 @@ long    _LikelihoodFunction::Bracket (long index, hyFloat& left, hyFloat& middle
 
     if (index < 0) {
         leftStep = middle;
+    } else {
+        if (middle - leftStep <= lowerBound) {
+            leftStep = MAX ((middle - lowerBound)*0.9, 1e-8);
+        }
     }
 
     if (saveM != middle) {
@@ -5330,7 +5334,7 @@ long    _LikelihoodFunction::Bracket (long index, hyFloat& left, hyFloat& middle
             leftStep = MIN (leftStep*0.125, middle-lowerBound);
             
             if ( leftStep<initialStep*.1 && index >= 0 || index < 0 && leftStep < STD_GRAD_STEP) {
-                if (!first || index >= 0 && current_log_l > -INFINITY) {
+                if (!first ||  index >= 0 && current_log_l > -INFINITY && leftStep <= 1e-8) {
                     if (go2Bound) {
                         middle = lowerBound;
                         
diff --git a/src/core/matrix.cpp b/src/core/matrix.cpp
index fe9e2619f..2692c7ac1 100644
--- a/src/core/matrix.cpp
+++ b/src/core/matrix.cpp
@@ -550,7 +550,7 @@ void    DuplicateMatrix (_Matrix* targetMatrix, _Matrix const* sourceMatrix) {
     
   targetMatrix->compressedIndex = nil;
   
-  if (sourceMatrix->theIndex) {
+  if (! sourceMatrix->is_dense()) {
     if (!(targetMatrix->theIndex = (long*)MatrixMemAllocate(sizeof(long) *sourceMatrix->lDim))) { // allocate element index storage
       HandleApplicationError ( kErrorStringMemoryFail );
     } else {
@@ -568,7 +568,7 @@ void    DuplicateMatrix (_Matrix* targetMatrix, _Matrix const* sourceMatrix) {
   targetMatrix->theData = nil;
   
   if (sourceMatrix->lDim) {
-    if (sourceMatrix->storageType==0)
+    if (sourceMatrix->is_polynomial())
       // matrix will store pointers to elements
     {
       if (targetMatrix->lDim) {
@@ -576,7 +576,7 @@ void    DuplicateMatrix (_Matrix* targetMatrix, _Matrix const* sourceMatrix) {
           HandleApplicationError ( kErrorStringMemoryFail );
         } else {
           memcpy ((void*)targetMatrix->theData,(void*)sourceMatrix->theData,sourceMatrix->lDim*sizeof(void*));
-          if (!sourceMatrix->theIndex) { // non-sparse matrix
+          if (sourceMatrix->is_dense()) { // non-sparse matrix
             for (long i=0; i<sourceMatrix->lDim; i++)
               if (sourceMatrix->GetMatrixObject(i)) {
                 (sourceMatrix->GetMatrixObject(i))->AddAReference();
@@ -591,12 +591,12 @@ void    DuplicateMatrix (_Matrix* targetMatrix, _Matrix const* sourceMatrix) {
           
         }
       }
-    } else if (sourceMatrix->storageType==2) {
+    } else if (sourceMatrix->is_expression_based()) {
       if (targetMatrix->lDim) {
         targetMatrix->theData = (hyFloat*)MatrixMemAllocate(sourceMatrix->lDim*sizeof(void*));
         _Formula ** theFormulas = (_Formula**)(sourceMatrix->theData), **newFormulas =
         (_Formula**)(targetMatrix->theData);
-        if (sourceMatrix->theIndex) {
+        if (sourceMatrix->is_dense() == false) {
           for (long i = 0; i<sourceMatrix->lDim; i++)
             if (sourceMatrix->IsNonEmpty(i)) {
               newFormulas[i] = (_Formula*)theFormulas[i]->makeDynamic();
@@ -653,10 +653,10 @@ _Matrix::_Matrix (long theHDim, long theVDim, bool sparse, bool allocateStorage)
 //_____________________________________________________________________________________________
 
 inline  bool    _Matrix::IsNonEmpty  (long logicalIndex) const {
-    if (theIndex) {
-        return theIndex[logicalIndex] != -1;
+    if (is_dense() == false) {
+        return theIndex [logicalIndex] != -1;
     }
-    if (storageType == _NUMERICAL_TYPE) {
+    if (is_numeric()) {
         return true;
     }
     return GetMatrixObject(logicalIndex)!=ZEROPOINTER;
@@ -707,7 +707,7 @@ bool        _Matrix::is_square_numeric(bool dense) const {
     if (storageType!=_NUMERICAL_TYPE || hDim != vDim || hDim==0L){
         HandleApplicationError ("Square numerical matrix required");
     }
-    if (dense && theIndex) {
+    if (dense && !is_dense()) {
         HandleApplicationError ("Square dense numerical matrix required");
     }
     return true;
@@ -1228,9 +1228,10 @@ HBLObjectRef   _Matrix::LUDecompose (void) const {
     }
     else {
       for (long i=0; i<lDim; i++) {
-        if (IsNonEmpty(i)) {
+        if (theIndex[i] != -1) {
           long cell_coord = theIndex[i];
-          result->Store(cell_coord/vDim,cell_coord%vDim,theData[i]);
+          long r = cell_coord / hDim;
+          result->Store(r,cell_coord-r*vDim,theData[i]);
         }
       }
     }
@@ -1491,37 +1492,77 @@ HBLObjectRef   _Matrix::MultByFreqs (long freqID, bool reuse_value_object) {
 
         if (theIndex) {
             _Matrix*    vm = (_Matrix*) value;
-            hyFloat *dp = vm ->theData;
-            hyFloat *tempDiags = (hyFloat*) alloca (sizeof(hyFloat) * hDim);
-            InitializeArray(tempDiags, hDim, 0.0);
-
-            if (freq_matrix) {
-                  for (long i=0; i<lDim; i++) {
-                      long p = theIndex[i];
-                      if (p != -1) {
-                          long h = p/vDim;
-                          p -= h*vDim;
-                          if (h!=p) {
-                              tempDiags[h] += (dp[i] *= freq_matrix->theData[p]);
+            hyFloat * __restrict dp = vm ->theData;
+
+            if (vm->compressedIndex) {
+                //vm->_validateCompressedStorage();
+                long from = 0L;
+                if (freq_matrix) {
+                    for (long r = 0; r < hDim; r++) {
+                        long diagEntry = -1;
+                        hyFloat diagAccumulator = 0.;
+                        for (long c = from; c < vm->compressedIndex[r]; c++) {
+                            long col_index = vm->compressedIndex[c+hDim];
+                            if (col_index != r) {
+                                dp[c] *= freq_matrix->theData[col_index];
+                                diagAccumulator -= dp[c];
+                            } else {
+                                diagEntry = c;
+                            }
+                        }
+                        from = vm->compressedIndex[r];
+                        dp[diagEntry] = diagAccumulator;
+                    }
+                } else {
+                    for (long r = 0; r < hDim; r++) {
+                        long diagEntry = -1;
+                        hyFloat diagAccumulator = 0.;
+                        for (long c = from; c < vm->compressedIndex[r]; c++) {
+                            //printf ("%ld\n", vm->theIndex[c]);
+                            if (vm->compressedIndex[c+hDim] != r) {
+                                diagAccumulator -= dp[c];
+                            } else {
+                                diagEntry = c;
+                            }
+                        }
+                        //printf ("%ld %ld %g\n", r, diagEntry, diagAccumulator);
+                        from = vm->compressedIndex[r];
+                        dp[diagEntry] = diagAccumulator;
+                    }
+                }
+                //vm->_validateCompressedStorage();
+            } else {
+                hyFloat *tempDiags = (hyFloat*) alloca (sizeof(hyFloat) * hDim);
+                InitializeArray(tempDiags, hDim, 0.0);
+
+                if (freq_matrix) {
+                      for (long i=0; i<lDim; i++) {
+                          long p = theIndex[i];
+                          if (p != -1) {
+                              long h = p / vDim;
+                                   p = p - h*vDim;
+                              if (h!=p) {
+                                  tempDiags[h] += (dp[i] *= freq_matrix->theData[p]);
+                              }
                           }
                       }
-                  }
-            }
-            else {
-                  for (long i=0; i<lDim; i++) {
-                      long p = theIndex[i];
-                      if (p != -1) {
-                          long h = p/vDim;
-                          p -= h*vDim;
-                          if (h!=p) {
-                              tempDiags[h] += dp[i];
+                }
+                else {
+                      for (long i=0; i<lDim; i++) {
+                          long p = theIndex[i];
+                          if (p != -1) {
+                              long h = p / vDim;
+                                   p = p - h*vDim;
+                              if (h!=p) {
+                                  tempDiags[h] += dp[i];
+                              }
                           }
                       }
-                  }
-            }
+                }
           
-            for (long j=0L; j<hDim; j++) {
-                vm->Store (j,j,-tempDiags[j]);
+                for (long j=0L; j<hDim; j++) {
+                    vm->Store (j,j,-tempDiags[j]);
+                }
             }
 
         } else {
@@ -1529,9 +1570,7 @@ HBLObjectRef   _Matrix::MultByFreqs (long freqID, bool reuse_value_object) {
 
             if (freq_matrix) {
                 if (freq_matrix->theIndex) {
-                    for (long i=0; i<lDim; i++) {
-                        theMatrix[i] *= (*freq_matrix)[i%vDim];
-                    }
+                    HandleApplicationError(_String("Sparse frequency matrices are not supported"));
                 } else {
                     for (unsigned long column=0UL; column<vDim; column++) {
                       const hyFloat freq_i = freq_matrix->theData[column];
@@ -1856,14 +1895,40 @@ bool    _Matrix::AmISparseFast (_Matrix& whereTo) {
     
     long * non_zero_index = (long*)alloca (threshold*sizeof(long));
     
-    for (unsigned long i=0UL; i<lDim ; i++) {
-          if (__builtin_expect(theData[i] != ZEROOBJECT, 0L)) {
-              non_zero_index[k++] = i;
-              if (k == threshold) {
-                  return false;
-              }
-          }
+#ifdef _SLKP_USE_AVX_INTRINSICS
+    __m256d zeros = _mm256_setzero_pd();
+    long lDimMOD4 = lDim >> 2 << 2;
+    for (long i = 0; i < lDimMOD4; i+=4) {
+         int res  = _mm256_movemask_pd(_mm256_cmp_pd (_mm256_loadu_pd (theData+i), zeros, _CMP_NEQ_OQ));
+         if (res) { // something is different
+            if (res & 1) { non_zero_index[k++] = i; if (k == threshold) break; };
+            if (res & 2) { non_zero_index[k++] = i+1; if (k == threshold) break; };
+            if (res & 4) { non_zero_index[k++] = i+2; if (k == threshold) break; };
+            if (res & 8) { non_zero_index[k++] = i+3; if (k == threshold) break; };
+        }
     }
+    
+    if (k < threshold)
+        for (long i = lDimMOD4; i < lDim; i++) {
+            if (theData[i] != 0.0) {
+                non_zero_index[k++] = i;
+                if (k == threshold) {
+                    return false;
+                }
+            }
+        }
+#else
+    for (long i = 0; i < lDim; i++) {
+        if (theData[i] != 0.0) {
+            non_zero_index[k++] = i;
+            if (k == threshold) {
+                return false;
+            }
+        }
+    }
+#endif
+   
+    
 
     if (k < threshold) {
         // we indeed are sparse enough
@@ -1888,14 +1953,15 @@ bool    _Matrix::AmISparseFast (_Matrix& whereTo) {
             whereTo.theIndex               = (long*)MatrixMemAllocate (whereTo.lDim*sizeof(long));
         }
         
-        hyFloat _hprestrict_ *          newData  =  canReuse ? whereTo.theData : (hyFloat*)MatrixMemAllocate (whereTo.lDim*sizeof(hyFloat));
+        hyFloat * _hprestrict_          newData  =  canReuse ? whereTo.theData : (hyFloat*)MatrixMemAllocate (whereTo.lDim*sizeof(hyFloat));
         
         if (whereTo.compressedIndex) {
             MatrixMemFree(whereTo.compressedIndex);
         }
         whereTo.compressedIndex = (long*) MatrixMemAllocate((whereTo.lDim + hDim) * sizeof (long));
         
-        long              currentRow = 0L;
+        long                    currentRow = 0L;
+        long   * __restrict     wci = (long*)whereTo.compressedIndex;
 
         for (long i=0L; i<k; i++) {
             //printf ("%ld %ld", i, theIndex[i]);
@@ -1905,10 +1971,10 @@ bool    _Matrix::AmISparseFast (_Matrix& whereTo) {
             long indexRow = entryIndex / vDim,
                  indexColumn = entryIndex - indexRow * vDim;
  
-            whereTo.compressedIndex[i + hDim] = indexColumn;
+            wci[i + hDim] = indexColumn;
             if (indexRow > currentRow) {
                 for (long l = currentRow; l < indexRow; l++) {
-                    whereTo.compressedIndex[l] = i;
+                    wci[l] = i;
                 }
                 currentRow = indexRow;
             }
@@ -2151,53 +2217,33 @@ bool    _Matrix::CheckIfSparseEnough(bool force, bool copy) {
 }
 
 //_____________________________________________________________________________________________
-bool    _Matrix::IncreaseStorage    (void)
-{
+bool    _Matrix::IncreaseStorage    (void) {
+    if (compressedIndex) {
+        HandleApplicationError("Internal error. Called _Matrix::IncreaseStorage on compressed index matrix");
+        return false;
+    }
     lDim += allocationBlock;
-
-    long* tempIndex, i;
-
-    if (!(tempIndex = (long*)MatrixMemAllocate(lDim*sizeof(long)))) {
-        HandleApplicationError ( kErrorStringMemoryFail );
-    } else {
-        memcpy (tempIndex, theIndex, (lDim-allocationBlock)*sizeof(long));
-        MatrixMemFree( theIndex);
-
-        for (i = lDim-1; i>=lDim-allocationBlock; i--) {
-            tempIndex [i] = -1;
-        }
-        theIndex = tempIndex;
+    theIndex = (long*)MemReallocate(theIndex,  lDim*sizeof(long));
+    for (long i = lDim-allocationBlock; i < lDim; i++) {
+        theIndex [i] = -1;
     }
 
-    if (storageType != 1)
+    if (!is_numeric()) {
         // pointers or formulas
-    {
         _MathObject** tempData;
         if (!(tempData = (_MathObject**) MatrixMemAllocate(sizeof( char)* lDim*sizeof(void*)))) {
             HandleApplicationError ( kErrorStringMemoryFail );
         } else {
-            memcpy (tempData, theData, (lDim-allocationBlock)*sizeof(void*));
-            MatrixMemFree (theData);
-            for (i = lDim-1; i>=lDim-allocationBlock; i--) {
-                tempData [i] = ZEROPOINTER;
+            theData = (hyFloat*) MemReallocate(theData, lDim*sizeof(void*));
+            for (long i = lDim-allocationBlock; i < lDim; i++) {
+                theData [i] = ZEROPOINTER;
             }
-            theData = (hyFloat*)tempData;
         }
-    } else
+    } else {
         //objects
-    {
-        hyFloat* tempData;
-        if (!(tempData =  (hyFloat*)MatrixMemAllocate(sizeof(hyFloat)* lDim))) {
-            HandleApplicationError ( kErrorStringMemoryFail );
-        } else {
-            for (i = lDim-1; i>=lDim-allocationBlock; i--) {
-                tempData [i] = ZEROOBJECT;
-            }
-            for (; i>=0; i--) {
-                tempData [i] = ((hyFloat*)theData) [i];
-            }
-            MatrixMemFree( theData);
-            theData = tempData;
+        theData = (hyFloat*) MemReallocate(theData, lDim*sizeof(hyFloat));
+        for (long i = lDim-allocationBlock; i < lDim; i++) {
+            theData [i] = ZEROOBJECT;
         }
     }
     return TRUE;
@@ -2209,10 +2255,10 @@ bool    _Matrix::IncreaseStorage    (void)
 
 void    _Matrix::Convert2Formulas (void)
 {
-    if (storageType == 1) {
-        storageType = 2;
+    if (is_numeric()) {
+        storageType = _FORMULA_TYPE;
         _Formula** tempData = (_Formula**)MatrixMemAllocate (sizeof(void*)*lDim);
-        if (!theIndex) {
+        if (is_dense()) {
             for (long i = 0; i<lDim; i++) {
                 tempData[i] = new _Formula (new _Constant (((hyFloat*)theData)[i]));
             }
@@ -2575,7 +2621,7 @@ _String*        _Matrix::BranchLengthExpression (_Matrix* baseFreqs, bool mbf) {
 
 //_____________________________________________________________________________________________
 void        _Matrix::MakeMeSimple (void) {
-    if (storageType == _FORMULA_TYPE) {
+    if (is_expression_based()) {
         long            stackLength = 0L;
 
         _SimpleList     newFormulas,
@@ -2588,6 +2634,10 @@ void        _Matrix::MakeMeSimple (void) {
         _SimpleList varListAux;
         _AVLList    varList (&varListAux);
       
+        
+        if (!is_dense()) {
+            CompressSparseMatrix(false, (hyFloat*)alloca (sizeof (hyFloat) * lDim));
+        }
         if (ProcessFormulas (stackLength,varList,newFormulas,references,flaStrings)) {
             storageType = _SIMPLE_FORMULA_TYPE;
 
@@ -2816,6 +2866,51 @@ void        _Matrix::FillInList (_List& fillMe, bool convert_numbers) const {
     }
 }
 
+//_____________________________________________________________________________________________
+bool   _Matrix::_validateCompressedStorage (void) const {
+    if (theIndex && compressedIndex) {
+        long from = 0L;
+        long last_index = -1L;
+        for (long r = 0; r < hDim; r++) {
+            if (compressedIndex[r] < from || compressedIndex[r] > lDim ) {
+                HandleApplicationError(_String ("Inconsistent compressedIndex row element count at " ) & r & " : " & from & " vs " & compressedIndex[r]);
+                return false;
+            }
+            for (long c = from; c < compressedIndex[r]; c++) {
+                long myIndex = theIndex[c];
+                if (myIndex <= last_index) {
+                    HandleApplicationError(_String ("Lack of sortedness in theIndex at " ) & c & " : " & myIndex & " vs " & last_index);
+                    return false;
+                }
+                if (myIndex < 0 || myIndex >= hDim * vDim) {
+                    HandleApplicationError(_String ("Out of bounds in theIndex at " ) & c & " : " & myIndex & ", lDim = " & lDim);
+                    return false;
+                }
+                last_index = myIndex;
+                if (c > from) {
+                    if (compressedIndex[c+hDim] <= compressedIndex [c-1+hDim]) {
+                        HandleApplicationError(_String ("Lack of sortedness in columns at " ) & c & " : " & compressedIndex[c+hDim] & " vs " & compressedIndex[c+hDim-1]);
+                        return false;
+                    }
+                }
+                if (myIndex / vDim != r || myIndex % hDim != compressedIndex[c+vDim]) {
+                    HandleApplicationError(_String ("Stored index does match row/column " ) & myIndex & "(" & c & ") : " & r & " , " & compressedIndex[c+hDim]);
+                    return false;
+                }
+            }
+            from = compressedIndex[r];
+        }
+        
+        if (compressedIndex[hDim-1] != lDim) {
+            HandleApplicationError(_String ("Incompatible compressedIndex[hDim-1] and lDim: " ) & compressedIndex[hDim-1] & " : " & lDim);
+            return false;
+        }
+        
+        return true;
+    }
+    return false;
+}
+
 //_____________________________________________________________________________________________
 HBLObjectRef   _Matrix::EvaluateSimple (_Matrix* existing_storage) {
 // evaluate the matrix  overwriting the old one
@@ -2848,46 +2943,152 @@ HBLObjectRef   _Matrix::EvaluateSimple (_Matrix* existing_storage) {
     }
 
 
-    for (long f = 0; f < cmd->formulasToEval.lLength; f++) {
+    for (long f = 0L; f < cmd->formulasToEval.lLength; f++) {
         cmd->formulaValues [f] = ((_Formula*)cmd->formulasToEval.list_data[f])->ComputeSimple(cmd->theStack, cmd->varValues);
     }
 
     long * fidx = cmd->formulaRefs;
 
     if (theIndex) {
-        if (result->lDim != lDim) {
-            result->lDim = lDim;
-            result->theIndex = (long*)MemReallocate((hyPointer)result->theIndex,sizeof(long)*lDim);
-            result->theData = (hyFloat*)MemReallocate ((hyPointer)result->theData,sizeof(hyFloat)*lDim);
-        }
+        
         result->bufferPerRow = bufferPerRow;
         result->overflowBuffer = overflowBuffer;
         result->allocationBlock = allocationBlock;
+        
+        long* diagIndices = nil;
 
+        if (compressedIndex) {
+            if (result->lDim < lDim + hDim) {
+                result->theIndex = (long*)MemReallocate((hyPointer)result->theIndex,sizeof(long)*(lDim + hDim));
+                result->theData = (hyFloat*)MemReallocate ((hyPointer)result->theData,sizeof(hyFloat)*(lDim+hDim));
+                
+            }
+            
+            result->lDim = lDim;
+            
+            if (result->compressedIndex) {
+                result->compressedIndex = (long*)MemReallocate ((hyPointer)result->compressedIndex,sizeof(long)*(lDim+hDim+hDim));
+            } else {
+                result->compressedIndex = (long*)MemAllocate (sizeof(long)*(lDim+hDim+hDim));
+            }
+            
+            if (hDim == vDim) {
+                
+                long elements_added                   = 0L;
+                long current_element_index_old_matrix = 0L;
+                long current_element_index_new_matrix = 0L;
+                long from = 0L;
+                auto copy_indices = [&] () -> void {
+                    result->theIndex[current_element_index_new_matrix] = theIndex[current_element_index_old_matrix];
+                    result->compressedIndex[current_element_index_new_matrix+hDim] = compressedIndex[hDim+current_element_index_old_matrix];
+                };
+                diagIndices = (long*)alloca (sizeof (long) * hDim);
+                auto inject_diagonal = [&] (long r) -> void {
+                    elements_added++;
+                    result->lDim ++;
+                    diagIndices [r] = current_element_index_new_matrix;
+                    result->theIndex[current_element_index_new_matrix] = r*vDim + r;
+                    result->theData[current_element_index_new_matrix] = 0.;
+                    result->compressedIndex[current_element_index_new_matrix+hDim] = r;
+                    current_element_index_new_matrix++;
+                };
+                /*if (!_validateCompressedStorage()) {
+                    HandleApplicationError("Error in compressed storage [before]");
+                }*/
+                for (long r = 0; r < hDim; r++) {
+                    diagIndices[r] = -1L;
+                    for (long c = from; c < compressedIndex[r]; c++, current_element_index_old_matrix++, current_element_index_new_matrix++) {
+                        if (compressedIndex[c + hDim] < r) { // column before diagonal; copy data
+                            result->theData[current_element_index_new_matrix] = cmd->formulaValues[fidx[current_element_index_old_matrix]];
+                            copy_indices();
+                        } else if (compressedIndex[c + hDim] > r) {
+                            if (diagIndices[r] == -1) { // no diagonal entry
+                                inject_diagonal(r);
+                            }
+                            result->theData[current_element_index_new_matrix] = cmd->formulaValues[fidx[current_element_index_old_matrix]];
+                            copy_indices();
+                        } else { // diagnoal entry
+                           copy_indices();
+                           diagIndices[r] = current_element_index_new_matrix;
+                        }
+                    }
+                    if (diagIndices[r] == -1) { // no diagonal entry
+                        inject_diagonal(r);
+                    }
+                    from = compressedIndex[r];
+                    result->compressedIndex[r] = from+elements_added;
+                    
+                }
+                /*if (!result->_validateCompressedStorage()) {
+                    HandleApplicationError("Error in compressed storage");
+                }*/
+            } else {
+                for (long i = 0; i<lDim; i++) {
+                    long idx = theIndex[i];
+                    result->theData[i] = cmd->formulaValues[fidx[i]];
+                    result->theIndex[i] = idx;
+                    result->compressedIndex[i] = compressedIndex[i];
+                }
+                for (long i = lDim; i<lDim+hDim; i++) {
+                    result->compressedIndex[i] = compressedIndex[i];
+                }
+            }
 
-        for (long i = 0; i<lDim; i++) {
-            long idx = theIndex[i];
-
-            if (idx != -1) {
-                result->theData[i] = cmd->formulaValues[fidx[i]];
+        } else {
+            
+            if (result->lDim != lDim) {
+                result->lDim = lDim;
+                result->theIndex = (long*)MemReallocate((hyPointer)result->theIndex,sizeof(long)*lDim);
+                result->theData = (hyFloat*)MemReallocate ((hyPointer)result->theData,sizeof(hyFloat)*lDim);
             }
+            
+            for (long i = 0; i<lDim; i++) {
+                long idx = theIndex[i];
+
+                if (idx != -1) {
+                    result->theData[i] = cmd->formulaValues[fidx[i]];
+                }
 
-            result->theIndex[i] = idx;
+                result->theIndex[i] = idx;
+            }
         }
 
 
         if (hDim==vDim) {
-            hyFloat* diagStorage = (hyFloat*)alloca (sizeof(hyFloat) * hDim);
-            InitializeArray(diagStorage, hDim, 0.0);
-            for (long i = 0; i<lDim; i++) {
-                long k = result->theIndex[i];
-                if (k!=-1) {
-                    diagStorage[k/hDim] += result->theData[i];
+            
+            if (result->compressedIndex) {
+                long from = 0L;
+                for (long r = 0; r < hDim; r++) {
+                    //printf ("%ld\n", diagIndices[r]);
+                    long di = diagIndices[r];
+                    for (long c = from; c < result->compressedIndex[r]; c++) {
+                        //printf ("%ld %g\n", c, result->theData[c]);
+                        if (c != di) {
+                            result->theData[di] -= result->theData[c];
+                        }
+                    }
+                    from = result->compressedIndex[r];
                 }
+                /*for (long r = 0; r < hDim; r++) {
+                    printf ("%ld %g\n", diagIndices[r], result->theData[diagIndices[r]]);
+                }
+                exit (0);*/
             }
-            for (long i = 0; i<hDim; i++) {
-                (*result)[i*hDim+i] = -diagStorage[i];
-            }
+            else {
+                hyFloat* diagStorage = (hyFloat*)alloca (sizeof(hyFloat) * hDim);
+                memset (diagStorage, 0, sizeof(hyFloat) * hDim);
+                for (long i = 0; i<lDim; i++) {
+                    long k = result->theIndex[i];
+                    if (k!=-1) {
+                        diagStorage[k/hDim] += result->theData[i];
+                    }
+                }
+                for (long i = 0; i<hDim; i++) {
+                     (*result)[i*hDim+i] = -diagStorage[i];
+                }
+           }
+            
+           
         }
     } else {
 
@@ -2997,7 +3198,7 @@ void    _Matrix::Clear (bool complete) {
 
 void    _Matrix::ZeroNumericMatrix (void) {
     if (is_numeric()) {
-        InitializeArray (theData, lDim, 0.0);
+        memset (theData, 0, sizeof (hyFloat) * lDim);
         if (!is_dense()) {
             InitializeArray (theIndex, lDim, -1L);
         }
@@ -3006,9 +3207,8 @@ void    _Matrix::ZeroNumericMatrix (void) {
 
 //_____________________________________________________________________________________________
 
-void    _Matrix::Resize (long newH)
-{
-    if (newH >= 0 && newH != hDim && storageType == 1 && theIndex == nil) {
+void    _Matrix::Resize (long newH) {
+    if (newH >= 0 && newH != hDim && is_numeric() && is_dense()) {
         hDim = newH;
         lDim = newH*vDim;
 
@@ -3151,13 +3351,24 @@ void    _Matrix::AddMatrix  (_Matrix& storage, _Matrix& secondArg, bool subtract
                        
             if (subtract) {
 #ifdef  _SLKP_USE_AVX_INTRINSICS
-            #define     CELL_OP(x) _mm256_storeu_pd (stData+x, _mm256_sub_pd (_mm256_loadu_pd (stData+x), _mm256_loadu_pd (argData+x)))
+            #define     CELL_OP1(x,y) __m256d y = _mm256_sub_pd (_mm256_loadu_pd (stData+x), _mm256_loadu_pd (argData+x))
+#define CELL_OP2(x,y) _mm256_storeu_pd (stData+x,y)
                 
-                for (long idx = 0; idx < upto; idx+=16) {
-                    CELL_OP (idx);
-                    CELL_OP (idx+4);
-                    CELL_OP (idx+8);
-                    CELL_OP (idx+12);
+        #pragma GCC unroll 4
+        #pragma clang loop vectorize(enable)
+        #pragma clang loop interleave(enable)
+        #pragma clang loop unroll(enable)
+        #pragma GCC ivdep
+        #pragma ivdep
+               for (long idx = 0; idx < upto; idx+=16) {
+                    CELL_OP1 (idx,r1);
+                    CELL_OP1 (idx+4,r2);
+                    CELL_OP1 (idx+8,r3);
+                    CELL_OP1 (idx+12,r4);
+                    CELL_OP2 (idx,r1);
+                    CELL_OP2 (idx+4,r2);
+                    CELL_OP2 (idx+8,r3);
+                    CELL_OP2 (idx+12,r4);
                 }
 #else
                 for (long idx = 0; idx < upto; idx+=4) {
@@ -3170,14 +3381,19 @@ void    _Matrix::AddMatrix  (_Matrix& storage, _Matrix& secondArg, bool subtract
             } else {
 #ifdef  _SLKP_USE_AVX_INTRINSICS
             #define     CELL_OP(x) _mm256_storeu_pd (stData+x, _mm256_add_pd (_mm256_loadu_pd (stData+x), _mm256_loadu_pd (argData+x)))
-                 
+
+
+            #pragma GCC unroll 4
+            #pragma clang loop vectorize(enable)
+            #pragma clang loop interleave(enable)
+            #pragma clang loop unroll(enable)
                  for (long idx = 0; idx < upto; idx+=16) {
                      CELL_OP (idx);
                      CELL_OP (idx+4);
                      CELL_OP (idx+8);
                      CELL_OP (idx+12);
                  }
-
+        
  #else
                 for (long idx = 0; idx < upto; idx+=4) {
                     stData[idx]+=argData[idx];
@@ -3493,12 +3709,14 @@ void    _Matrix::Multiply  (_Matrix& storage, _Matrix const& secondArg) const
                 hyFloat  * dest = storage.theData;
                 
 #if defined _SLKP_USE_AVX_INTRINSICS
+                #define DO_GROUP_OP0(X,Y,k) Y = _mm256_loadu_pd(secondArg.theData + col_offset + k); X = _mm256_mul_pd (A4,Y);
                 #ifdef _SLKP_USE_FMA3_INTRINSICS
                     #define DO_GROUP_OP(X,Y,k) X = _mm256_loadu_pd(dest + row_offset + k); Y = _mm256_loadu_pd(secondArg.theData + col_offset + k); _mm256_storeu_pd (dest + row_offset + k, _mm256_fmadd_pd (A4,Y,X));
                     #define DO_GROUP_OP1(X,Y,k) X = _mm256_loadu_pd(dest + row_offset + k); Y = _mm256_loadu_pd(secondArg.theData + col_offset + k); X = _mm256_fmadd_pd (A4,Y,X);
                     #define DO_GROUP_OP2(X,k) _mm256_storeu_pd (dest + row_offset + k,X);
                 #else
                     #define DO_GROUP_OP(X,Y,k) X = _mm256_loadu_pd(dest + row_offset + k);  Y = _mm256_loadu_pd(secondArg.theData + col_offset + k); _mm256_storeu_pd (dest + row_offset + k, _mm256_add_pd (X, _mm256_mul_pd(A4, Y)));
+
                     #define DO_GROUP_OP1(X,Y,k) X = _mm256_loadu_pd(dest + row_offset + k); Y = _mm256_loadu_pd(secondArg.theData + col_offset + k); X = _mm256_add_pd (X, _mm256_mul_pd(A4, Y));
                     #define DO_GROUP_OP2(X,k) _mm256_storeu_pd (dest + row_offset + k,X);
                #endif
@@ -3513,7 +3731,6 @@ void    _Matrix::Multiply  (_Matrix& storage, _Matrix const& secondArg) const
             
               
               if (dimm4 == vDim) {
-                  memset (dest, 0, lDim * sizeof (hyFloat));
                   #if defined  _SLKP_USE_AVX_INTRINSICS
                       long ti = 0L,
                            row_offset = 0L;
@@ -3521,10 +3738,23 @@ void    _Matrix::Multiply  (_Matrix& storage, _Matrix const& secondArg) const
                       if (vDim == 20UL) {
                           for (long r = 0; r < 20; r++, row_offset += 20) {
                               long col_offset = 0L;
-                              for (long c = 0; c < 20L; c++, ti++, col_offset += 20L) {
-                                  __m256d A4 = _mm256_set1_pd(theData[ti]);
+                              __m256d A4 = _mm256_set1_pd(theData[ti]);
+                              __m256d D4, B4, D4_1, B4_1, D4_2, B4_2, D4_3, B4_3, D4_4, B4_4;
+                              DO_GROUP_OP0 (D4, B4, 0);
+                              DO_GROUP_OP0 (D4_1, B4_1, 4);
+                              DO_GROUP_OP0 (D4_2, B4_2, 8);
+                              DO_GROUP_OP0 (D4_3, B4_3, 12);
+                              DO_GROUP_OP0 (D4_4, B4_4, 16);
+                              DO_GROUP_OP2 (D4, 0);
+                              DO_GROUP_OP2 (D4_1, 4);
+                              DO_GROUP_OP2 (D4_2, 8);
+                              DO_GROUP_OP2 (D4_3, 12);
+                              DO_GROUP_OP2 (D4_4, 16);
+                              ti++;
+                              col_offset = 20L;
+                              for (long c = 1; c < 20L; c++, ti++, col_offset += 20L) {
+                                  A4 = _mm256_set1_pd(theData[ti]);
                                   //for (long k = 0; k < 20L; k+=4) {
-                                      __m256d D4, B4, D4_1, B4_1, D4_2, B4_2, D4_3, B4_3, D4_4, B4_4;
                                       DO_GROUP_OP1 (D4, B4, 0);
                                       DO_GROUP_OP1 (D4_1, B4_1, 4);
                                       DO_GROUP_OP1 (D4_2, B4_2, 8);
@@ -3541,8 +3771,16 @@ void    _Matrix::Multiply  (_Matrix& storage, _Matrix const& secondArg) const
                       } else if (vDim == 60UL) {
                           for (long r = 0; r < 60; r++, row_offset += 60) {
                               long col_offset = 0L;
-                              for (long c = 0; c < 60L; c++, ti++, col_offset += 60L) {
-                                  __m256d A4 = _mm256_set1_pd(theData[ti]);
+                              __m256d A4 = _mm256_set1_pd(theData[ti]);
+                              for (long k = 0; k < 60L; k+=4) {
+                                  __m256d D4, B4;
+                                  DO_GROUP_OP0 (D4, B4, k);
+                                  DO_GROUP_OP2 (D4, k);
+                              }
+                              col_offset = 60L;
+                              ti++;
+                              for (long c = 1; c < 60L; c++, ti++, col_offset += 60L) {
+                                  A4 = _mm256_set1_pd(theData[ti]);
                                   for (long k = 0; k < 60L; k+=4) {
                                       __m256d D4, B4;
                                       DO_GROUP_OP (D4, B4, k);
@@ -3552,8 +3790,16 @@ void    _Matrix::Multiply  (_Matrix& storage, _Matrix const& secondArg) const
                       } else
                           for (long r = 0; r < vDim; r++, row_offset += vDim) {
                               long col_offset = 0L;
-                              for (long c = 0; c < vDim; c++, ti++, col_offset += vDim) {
-                                  __m256d A4 = _mm256_set1_pd(theData[ti]);
+                              __m256d A4 = _mm256_set1_pd(theData[ti]);
+                              for (long k = 0; k < vDim; k+=4) {
+                                  __m256d D4, B4;
+                                  DO_GROUP_OP0 (D4, B4, k);
+                                  DO_GROUP_OP2 (D4, k);
+                              }
+                              col_offset = vDim;
+                              ti++;
+                              for (long c = 1; c < vDim; c++, ti++, col_offset += vDim) {
+                                 A4 = _mm256_set1_pd(theData[ti]);
                                  #pragma GCC unroll 4
                                  #pragma clang loop vectorize(enable)
                                  #pragma clang loop interleave(enable)
@@ -3566,7 +3812,8 @@ void    _Matrix::Multiply  (_Matrix& storage, _Matrix const& secondArg) const
                           }
                       return;
                   #elif _SLKP_USE_SSE_INTRINSICS
-                      long ti = 0L,
+                    memset (dest, 0, lDim * sizeof (hyFloat));
+                     long ti = 0L,
                            row_offset = 0L;
                   
                       if (vDim == 20UL) {
@@ -3617,7 +3864,8 @@ void    _Matrix::Multiply  (_Matrix& storage, _Matrix const& secondArg) const
                           }
                       return;
                 #endif
-                  
+                memset (dest, 0, lDim * sizeof (hyFloat));
+
                 for (unsigned long c = 0UL; c < secondArg.vDim; c ++) {
 
                   /*
@@ -4481,7 +4729,42 @@ hyFloat  _Matrix::MaxElement  (char runMode, long* indexStore) const {
             }
             return max;
         } else {
-            for (long i = 0; i<lDim; i++) {
+            
+            if (doAbsValue) {
+                if (doMaxElement) {
+                    for (long i = 0; i<lDim; i++) {
+                        hyFloat t = fabs(theData[i]);
+                        if (t > max) {
+                            max = t;
+                            if (indexStore) {
+                                *indexStore = i;
+                            }
+                        }
+                    }
+                } else {
+                    for (long i = 0; i<lDim; i++) {
+                        max += fabs (theData[i]);
+                    }
+                }
+            } else {
+                if (doMaxElement) {
+                    for (long i = 0; i<lDim; i++) {
+                        hyFloat t = theData[i];
+                        if (t > max) {
+                            max = t;
+                            if (indexStore) {
+                                *indexStore = i;
+                            }
+                        }
+                    }
+                } else {
+                    for (long i = 0; i<lDim; i++) {
+                        max += theData[i];
+                    }
+                }
+            }
+            
+            /*for (long i = 0; i<lDim; i++) {
                 temp = theData[i];
                 if (doAbsValue && temp<0.0) {
                     temp = -temp;
@@ -4497,7 +4780,7 @@ hyFloat  _Matrix::MaxElement  (char runMode, long* indexStore) const {
                 } else {
                     max += temp;
                 }
-            }
+            }*/
             return max;
         }
     }
@@ -4510,13 +4793,11 @@ hyFloat  _Matrix::MaxElement  (char runMode, long* indexStore) const {
 
 //_____________________________________________________________________________________________
 
-void    _Matrix::RowAndColumnMax  (hyFloat& r, hyFloat &c, hyFloat * cache)
+void    _Matrix::RowAndColumnMax  (hyFloat& r, hyFloat &c, hyFloat * cache) {
 
 // returns the maximum row sum / column sum
 // the cache must be big enough to hold hDim + vDim
 // leave as nil to allocate cache run time
-
-{
     r = c = 10.;
 
     if (is_numeric()) { // numeric matrix
@@ -4531,23 +4812,39 @@ void    _Matrix::RowAndColumnMax  (hyFloat& r, hyFloat &c, hyFloat * cache)
         hyFloat * rowMax = maxScratch,
                      * colMax = maxScratch + hDim;
 
-        if (theIndex)
-            // sparse matrix
-            for (long i = 0; i<lDim; i++) {
-                long k = theIndex[i];
-                if  (k!=-1) {
-                    hyFloat temp = theData[i];
+        if (theIndex) {
+            if (compressedIndex) {
+                long from = 0L;
+                for (long r = 0; r < hDim; r++) {
+                    for (long c = from; c < compressedIndex[r]; c++) {
+                        hyFloat temp = theData[c];
+                        if (temp > 0.0) {
+                            rowMax[r] += temp;
+                            colMax[compressedIndex[c+hDim]] += temp;
+                        } else {
+                            rowMax[r] -= temp;
+                            colMax[compressedIndex[c+hDim]] -= temp;
+                        }
+                    }
+                    from = compressedIndex[r];
+                }
+            } else {
+                for (long i = 0; i<lDim; i++) {
+                    long k = theIndex[i];
+                    if  (k!=-1) {
+                        hyFloat temp = theData[i];
 
-                    if (temp<0.0) {
-                        rowMax[k/vDim] -= temp;
-                        colMax[k%vDim] -= temp;
-                    } else {
-                        rowMax[k/vDim] += temp;
-                        colMax[k%vDim] += temp;
+                        if (temp>0.0) {
+                            rowMax[k/vDim] += temp;
+                            colMax[k%vDim] += temp;
+                        } else {
+                            rowMax[k/vDim] -= temp;
+                            colMax[k%vDim] -= temp;
+                        }
                     }
                 }
             }
-        else
+        } else
             // dense matrix
             for (long i = 0, k=0; i<hDim; i++) {
                 for (long j=0; j<vDim; j++, k++) {
@@ -4585,7 +4882,7 @@ bool    _Matrix::IsMaxElement  (hyFloat bench) {
         if (!theIndex || compressedIndex) {
             for (long i = 0; i<lDim; i++) {
                 hyFloat t = theData[i];
-                if ( t<mBench || t>bench ) {
+                if ( t>bench || t<mBench ) {
                     return true;
                 }
             }
@@ -4593,7 +4890,7 @@ bool    _Matrix::IsMaxElement  (hyFloat bench) {
             for (long i = 0; i<lDim; i++) {
                 if (theIndex [i] >= 0) {
                     hyFloat t = theData[i];
-                    if ( t<mBench || t>bench ) {
+                    if ( t>bench || t<mBench ) {
                         return true;
                     }
                 }
@@ -4841,6 +5138,73 @@ void    _Matrix::CompressSparseMatrix (bool transpose, hyFloat * stash) {
     
     
     if (theIndex) {
+        
+        if (compressedIndex && transpose) {
+            long ai = 0;
+            long from = 0L;
+            memcpy (stash, theData, sizeof (hyFloat)*lDim);
+            long * by_column_counts = (long*)alloca ((hDim + lDim) * sizeof (long));
+            
+            memset (by_column_counts, 0, hDim * sizeof (long));
+            
+            // pass 1 : count records by column
+            long * stored_by_col = (long*)alloca (hDim * sizeof (long));
+            memset (stored_by_col, 0, hDim * sizeof (long));
+
+            for (long r = 0L; r < hDim; r++) {
+                for (long c = from; c < compressedIndex[r]; c++) {
+                    by_column_counts[compressedIndex[hDim+c]] ++;
+                }
+                from = compressedIndex[r];
+            }
+            
+            
+            for (long r = 1L; r < hDim; r++) {
+                by_column_counts[r] += by_column_counts[r-1];
+                //printf ("%ld %ld (%ld)\n", r, by_column_counts[r], compressedIndex[r]);
+            }
+            
+            // pass 2 : invert indices
+            from = 0L;
+            
+            
+            for (long r = 0L; r < hDim; r++) {
+                for (long c = from; c < compressedIndex[r]; c++) {
+                    long my_col = compressedIndex[c+hDim];
+                    // (i,j) => (j,i)
+                    long new_index = my_col * hDim + r;
+                    long col_offset = my_col ? by_column_counts[my_col-1] : 0;
+                    long new_compressed_index = stored_by_col[my_col] + col_offset;
+                    //printf ("%ld => %ld (%ld, %ld)\n", c, new_compressed_index,r*hDim + my_col, new_index);
+                    theIndex[new_compressed_index] = new_index;
+                    theData[new_compressed_index] = stash[c];
+                    by_column_counts[hDim + new_compressed_index] = r;
+                    stored_by_col[my_col]++;
+                }
+                from = compressedIndex[r];
+            }
+            
+            
+            
+            memcpy (compressedIndex, by_column_counts, sizeof (long)*(lDim+hDim));
+            //_validateCompressedStorage();
+            
+            /*for (long r = 0; r < hDim; r++) {
+                long trow = r;
+                for (long c = from; c < compressedIndex[r]; c++) {
+                    long tcol = compressedIndex[hDim+c];
+                    long transposed_index = tcol * vDim + trow;
+                    stash[ai] = theData[c];
+                    //sortedIndex.list_data[ai] = transposed_index;
+                    //sortedIndex3.list_data[ai] = transposed_index;
+                    ai++;
+                    //if (max < transposed_index) max = transposed_index;
+                }
+                from = compressedIndex[r];
+            }*/
+            return;
+        }
+        
         _SimpleList sortedIndex  ((unsigned long)lDim, (long*)alloca (lDim * sizeof (long))),
                     sortedIndex3 ((unsigned long)lDim, (long*)alloca (lDim * sizeof (long))),
                     sortedIndex2;
@@ -4855,44 +5219,27 @@ void    _Matrix::CompressSparseMatrix (bool transpose, hyFloat * stash) {
 
 
         if (transpose) {
+            long ai = 0L;
+            
             for (long i2=0; i2<lDim; i2++) {
                 long k = theIndex[i2];
                 if  (__builtin_expect (k!=-1,1)) {
-                    /*
-                    long r = k/vDim,
-                         c = k%vDim,
-                         r2 = c / blockChunk * blockShift + r / blockChunk,
-                         r3 = r2 * lDim + r * vDim + c;
-
-                    sortedIndex  << (c*vDim + r);
-                    sortedIndex3 << r3;
-                    stash[sortedIndex.lLength-1] = theData[i2];
-                    if (r3 > max) {
-                        max = r3;
-                    }*/
-                    long transposed_index = (k%vDim) * vDim + k/vDim;
-                    stash[sortedIndex.lLength] = theData[i2];
-                    sortedIndex  << transposed_index;
-                    sortedIndex3 << transposed_index;
+                    long trow = k/vDim, tcol = k - trow*vDim;
+                    long transposed_index = tcol * vDim + trow;
+                    stash[ai] = theData[i2];
+                    sortedIndex.list_data[ai] = transposed_index;
+                    sortedIndex3.list_data[ai] = transposed_index;
+                    ai++;
                     if (max < transposed_index) max = transposed_index;
-                    
                 }
             }
+            
+            sortedIndex3.lLength = ai;
+            sortedIndex.lLength = ai;
         } else {
             for (long i2=0; i2<lDim; i2++) {
                 long k = theIndex[i2];
                 if  (__builtin_expect (k!=-1,1)) {
-                    /*long r = k%vDim,
-                         c = k/vDim,
-                         r2 = c / blockChunk * blockShift + r / blockChunk,
-                         r3 = r2 * lDim + r * vDim + c;
-
-                    sortedIndex  << (c*vDim + r);
-                    sortedIndex3 << r3;
-                    stash[sortedIndex.lLength-1] = theData[i2];
-                    if (r3 > max) {
-                        max = r3;
-                    }*/
                     stash[sortedIndex.lLength] = theData[i2];
                     sortedIndex  << k;
                     sortedIndex3 << k;
@@ -4957,32 +5304,6 @@ void    _Matrix::CompressSparseMatrix (bool transpose, hyFloat * stash) {
         
         for (long l = currentRow; l < firstDim; l++)
             compressedIndex[l] = lDim;
-        //compressedIndex[firstDim-1] = lDim;
-        //printf (">[%ld] %ld\n", firstDim-1, compressedIndex[firstDim-1]);
-        //exit (1);
-        
-        
-        // this code checks conversion consistency
-         
-        /*long ei = 0;
-        long from = 0;
-        for (long i=0; i<hDim; i++) {
-            for (long k = from; k < compressedIndex[i]; k++, ei++) {
-                long my_row = theIndex[ei] / vDim,
-                     my_column = theIndex[ei] % vDim;
-                
-                if (my_row != i || my_column != compressedIndex[hDim+k]) {
-                    printf ("BARF\n");
-                }
-            }
-            from = compressedIndex[i];
-        }
-        
-        if (ei != lDim) {
-            printf ("BARF\n");
-        }*/
-        
-
     }
 }
 
@@ -5202,11 +5523,27 @@ _Matrix*    _Matrix::Exponentiate (hyFloat scale_to, bool check_transition, _Mat
         
         if (theIndex) {
             // transpose back
-            for (i=0; i<lDim; i++) {
-                long k = theIndex[i];
-                if  (k!=-1) {
-                    long div = k / vDim;
-                    theIndex[i] = (k - div * vDim)*vDim + div;
+            if (compressedIndex) {
+                long from = 0L, i = 0;
+                for (long r = 0; r < hDim; r++) {
+                    #pragma GCC unroll 4
+                    #pragma clang loop vectorize(enable)
+                    #pragma clang loop interleave(enable)
+                    #pragma clang loop unroll(enable)
+                    for (long c = from; c < compressedIndex[r]; c++, i++) {
+                        theIndex[i] = compressedIndex[c+hDim] * vDim + r;
+                    }
+                    from = compressedIndex[r];
+                }
+                MatrixMemFree(compressedIndex);
+                compressedIndex = nil;
+            } else {
+                for (i=0; i<lDim; i++) {
+                    long k = theIndex[i];
+                    if  (k!=-1) {
+                        long div = k / vDim;
+                        theIndex[i] = (k - div * vDim)*vDim + div;
+                    }
                 }
             }
             result->Transpose();
@@ -5323,12 +5660,41 @@ long    _Matrix::Hash  (long i, long j) const {
         return i*vDim+j;
     }
     // ordinary matrix
+    
+    if (compressedIndex) {
+        for (long c = i ? compressedIndex[i-1] : 0; c < compressedIndex[i]; c++) {
+            if (compressedIndex[hDim+c] == j) return c;
+        }
+        return -1;
+    }
 
-    long elementIndex = i*vDim+j, k, l, m=i*bufferPerRow,n,p;
+    long elementIndex   = i*vDim+j,
+         m              = i*bufferPerRow;
+    
+    for (long allocation_block_index = 0L; allocation_block_index < lDim; allocation_block_index += allocationBlock, m += allocationBlock) {
+        // look within the row for this allocation block
+        for (long l = m; l < m + bufferPerRow; l++) {
+            long try_me = theIndex[l];
+            if (try_me != elementIndex) {
+                if (try_me == -1) return -l - 2;
+            } else return l;
+        }
+        // if not found, look in the overflow are for this block
+        long upper_bound = MIN (lDim, allocation_block_index + allocationBlock);
+        for (long n = upper_bound - overflowBuffer; n < upper_bound; n++) {
+            long try_me = theIndex[n];
+            if (try_me != elementIndex) {
+                if (try_me == -1) return -n - 2;
+            } else return n;
+        }
+    }
+    return -1;
+    
 
-    for (k = 0; k<lDim/allocationBlock; k++,m+=allocationBlock) {
-        for (l=m; l<m+bufferPerRow; l++) {
-            p = theIndex[l];
+    
+    /*for (long blockIndex = 0; blockIndex<lDim/allocationBlock; blockIndex++,m+=allocationBlock) {
+        for (long l=m; l<m+bufferPerRow; l++) {
+            long p = theIndex[l];
             if (p!=elementIndex) {
                 if (p==-1) {
                     return -l-2;
@@ -5337,9 +5703,9 @@ long    _Matrix::Hash  (long i, long j) const {
                 return l;
             }
         }
-        n = (k+1)*allocationBlock-1;
-        for (l = n; l>n-overflowBuffer; l--) {
-            p = theIndex[l];
+        long n = (blockIndex+1)*allocationBlock-1;
+        for (long l = n; l>n-overflowBuffer; l--) {
+            long p = theIndex[l];
             if (p!=elementIndex) {
                 if (p==-1) {
                     return -l-2;
@@ -5348,8 +5714,8 @@ long    _Matrix::Hash  (long i, long j) const {
                 return l;
             }
         }
-    }
-    return -1;
+    }*/
+    //return -1;
 }
 
 //_____________________________________________________________________________________________
diff --git a/src/core/operation.cpp b/src/core/operation.cpp
index 6bfdf7902..b4cf8c880 100644
--- a/src/core/operation.cpp
+++ b/src/core/operation.cpp
@@ -662,10 +662,10 @@ bool        _Operation::ExecutePolynomial (_Stack& theScrap, _VariableContainer*
         } else {
             return false;
         }
-    }
-
-    if (theData>-1) {
-        p=  new _Polynomial(*LocateVar(theData>-1?theData:-theData-2));
+    } else {
+        if (theData>-1) {
+            p=  new _Polynomial(*LocateVar(theData>-1?theData:-theData-2));
+        }
     }
 
     if (p) {
diff --git a/src/core/tree.cpp b/src/core/tree.cpp
index 759c55b33..35e8feadd 100644
--- a/src/core/tree.cpp
+++ b/src/core/tree.cpp
@@ -3050,6 +3050,11 @@ hyFloat          _TheTree::ComputeLLWithBranchCache (
             storageVec [direct_index] = accumulator;
         } else {
             if (accumulator <= 0.0) {
+                //fprintf (stderr, "ZERO TERM AT SITE %ld (direct %ld) EVAL %ld\n",siteID,direct_index, likeFuncEvalCallCount);
+                /*for (long s = 0; s < theFilter->NumberSpecies(); s++) {
+                    fprintf (stderr, "%s", theFilter->RetrieveState(direct_index, s).get_str());
+                }
+                fprintf (stderr, "\n");*/
                 throw (1L+direct_index);
             }
             
@@ -3061,7 +3066,7 @@ hyFloat          _TheTree::ComputeLLWithBranchCache (
                 term = log(accumulator) - correction;
             }
             
-            /*if (likeFuncEvalCallCount == 11035) {
+            /*if (likeFuncEvalCallCount == 15098) {
                 fprintf (stderr, "CACHE, %ld, %ld, %20.15lg, %20.15lg, %20.15lg\n", likeFuncEvalCallCount, siteID, accumulator, correction, term);
             }*/
             
@@ -3093,6 +3098,12 @@ hyFloat          _TheTree::ComputeLLWithBranchCache (
     
     
     // cases by alphabet dimension
+    /*if (likeFuncEvalCallCount == 15098) {
+        fprintf (stderr, "\nBRANCH ID %ld (%s)\n", brID, givenTreeNode->GetName()->get_str());
+        for (long e = 0; e < 16; e++) {
+            fprintf (stderr, "%ld => %lg\n", transitionMatrix[e]);
+        }
+    }*/
     
     try {
         switch (alphabetDimension) {
diff --git a/src/core/tree_evaluator.cpp b/src/core/tree_evaluator.cpp
index b786a28d7..0b1d316d8 100644
--- a/src/core/tree_evaluator.cpp
+++ b/src/core/tree_evaluator.cpp
@@ -45,6 +45,8 @@
 #include "global_things.h"
 #include "likefunc.h"
 
+extern  long likeFuncEvalCallCount;
+
 
 using namespace hy_global;
 using namespace hy_env;
@@ -144,11 +146,24 @@ template<long D> inline bool __ll_handle_conditional_array_initialization ( long
         }
         if (__builtin_expect(siteState >= 0L,1)) {
             // a single character state; sweep down the appropriate column
+            /*if (likeFuncEvalCallCount == 15098 && nodeCode == 3706 && siteID == 91) {
+                fprintf (stderr, "\nSITE CHECK: ID %ld, STATE %ld\n", nodeCode, siteState);
+                for (long e = 0; e < 4; e++) {
+                    fprintf (stderr, "%ld => %lg\n", e, parentConditionals[e]);
+                }
+            }*/
+            
             #pragma unroll(4)
             #pragma GCC unroll 4
             for (long k = 0L; k < D; k++) {
                 parentConditionals[k] *= tMatrix[siteState+D*k];
             }
+            /*if (likeFuncEvalCallCount == 15098 && nodeCode == 3706 && siteID == 91) {
+                fprintf (stderr, "\nSITE CHECK: ID %ld, STATE %ld\n", nodeCode, siteState);
+                for (long e = 0; e < 4; e++) {
+                    fprintf (stderr, "%ld => %lg\n", e, parentConditionals[e]);
+                }
+            }*/
             return true;
         } else {
             childVector = lNodeResolutions->theData + (-siteState-1) * D;
@@ -471,11 +486,18 @@ template<long D, bool ADJUST> inline void __ll_loop_handle_scaling (hyFloat& sum
         
         hyFloat scaler = _computeBoostScaler(scalingAdjustments [parentCode*siteCount + siteID] * _lfScalerUpwards, sum, didScale);
         
+        
+        /*if (likeFuncEvalCallCount == 15098 && siteID == 91) {
+            fprintf (stderr, "UP %ld (%ld) %lg\n", didScale, parentCode, scaler);
+        }*/
         if (didScale) {
             #pragma unroll(4)
             #pragma GCC unroll 4
             for (long c = 0; c < D; c++) {
                 parentConditionals [c] *= scaler;
+                /*if (likeFuncEvalCallCount == 15098 && siteID == 91) {
+                    fprintf (stderr, "%ld=>%g\n", c, parentConditionals [c]);
+                }*/
             }
             
             if (siteFrequency == 1L) {
@@ -492,13 +514,19 @@ template<long D, bool ADJUST> inline void __ll_loop_handle_scaling (hyFloat& sum
             if (sum < HUGE_VAL) { // no point scaling an infinity
                 
                 hyFloat scaler = _computeReductionScaler (scalingAdjustments [parentCode*siteCount + siteID] * _lfScalingFactorThreshold, sum, didScale);
+                /*if (likeFuncEvalCallCount == 15098 && siteID == 91) {
+                    fprintf (stderr, "DOWN %ld (%ld) %lg\n", didScale, parentCode, scaler);
+                }*/
                 
                 if (didScale) {
                     #pragma unroll(4)
                     #pragma GCC unroll 4
                     for (long c = 0; c < D; c++) {
-                         parentConditionals [c] *= scaler;
-                     }
+                        parentConditionals [c] *= scaler;
+                        /*if (likeFuncEvalCallCount == 15098 && siteID == 91) {
+                            fprintf (stderr, "%ld=>%g\n", c, parentConditionals [c]);
+                        }*/
+                    }
                     
                     if (siteFrequency == 1L) {
                         localScalerChange += didScale;
@@ -858,7 +886,16 @@ hyFloat      _TheTree::ComputeTreeBlockByBranch  (                   _SimpleList
         
         hyFloat  const * _hprestrict_ transitionMatrix = currentTreeNode->GetCompExp(catID)->theData;
         
-        
+        /*
+         if (likeFuncEvalCallCount == 15098 && parentCode == 3688) {
+            //if (currentTreeNode->GetName()->EndsWith("mt811400_SARS2_orf1ab_usa__4")) {
+            fprintf (stderr, "\nBRANCH ID %ld (%ld = parent ID, parent name = %s) (%s)\n", nodeCode, parentCode, ((_CalcNode*)flatTree (parentCode))->GetName()->get_str(), currentTreeNode->GetName()->get_str());
+                for (long e = 0; e < 16; e++) {
+                    fprintf (stderr, "%ld => %lg\n", transitionMatrix[e]);
+                }
+            //}
+        }
+         */
         
         long currentTCCIndex,currentTCCBit,parentTCCIIndex,parentTCCIBit;
         __ll_handle_tcc_init (tcc, isLeaf, siteCount, siteFrom, nodeCode, parentCode, parentTCCIBit, parentTCCIIndex, currentTCCBit, currentTCCIndex);
@@ -885,7 +922,8 @@ hyFloat      _TheTree::ComputeTreeBlockByBranch  (                   _SimpleList
                                                                       lNodeFlags, isLeaf, nodeCode, setBranch, flatTree.lLength, siteID, siteFrom, siteCount, siteOrdering, parentConditionals, tMatrix, lNodeResolutions, childVector, tcc, currentTCCBit, currentTCCIndex, lastUpdatedSite, setBranchTo)) {
                     continue;
                 }
-
+                
+                
                      
                 #ifdef _SLKP_USE_AVX_INTRINSICS
                     _handle4x4_pruning_case (childVector, tMatrix, parentConditionals, tmatrix_transpose);
@@ -1245,7 +1283,9 @@ hyFloat      _TheTree::ComputeTreeBlockByBranch  (                   _SimpleList
                 accumulator += rootConditionals[rootIndex] * theProbs[p];
             }
         }
-           
+        /*if (likeFuncEvalCallCount == 15098 && siteID == 91) {
+            fprintf (stderr, "\nREGULAR COMPUTE %lg (%ld) (%lg %lg %lg %lg)\n", accumulator, setBranch, rootConditionals[rootIndex-4],rootConditionals[rootIndex-3],rootConditionals[rootIndex-2],rootConditionals[rootIndex-1]);
+        }*/
         if (storageVec) {
             storageVec [siteOrdering.list_data[siteID]] = accumulator;
         } else {
@@ -1297,7 +1337,6 @@ template<long D> inline bool __lcache_loop_preface (bool isLeaf, long* __restric
             for (long k = 0L; k < D; k++, target_index+=D) {
                 parentConditionals[k]   *= tMatrix[target_index];
             }
-            
             return true;
         } else {
             childVector = lNodeResolutions->theData + (-siteState-1) * D;
@@ -1563,6 +1602,16 @@ void            _TheTree::ComputeBranchCache    (
     taggedNodes.Populate (flatTree.lLength, 0, 0);
     rootPath.Flip ();
     
+    /*if (likeFuncEvalCallCount == 15098) {
+        fprintf (stderr, "\nCaching branch %ld\n", brID);
+        rootPath.Each ([&](long n, unsigned long i) -> void {
+            _CalcNode * current_node = (_CalcNode*) (n < flatLeaves.lLength ?  flatCLeaves (n):
+                                                                              flatTree    (n-flatLeaves.lLength));
+            fprintf (stderr, "[%d: %d] %s\n", i, n,current_node->GetName()->get_str());
+        });
+        ObjectToConsole(&nodesToProcess);NLToConsole();
+    }*/
+    
     long const node_count = nodesToProcess.lLength + rootPath.lLength - 2L;
     
     for  (long nodeID = 0; nodeID < node_count; nodeID++) {
@@ -1585,9 +1634,8 @@ void            _TheTree::ComputeBranchCache    (
         }
         
         hyFloat * parentConditionals = iNodeCache +            (siteFrom + parentCode  * siteCount) * alphabetDimension;
-        if (taggedNodes.list_data[parentCode] == 0L)
+        if (taggedNodes.list_data[parentCode] == 0L) {
             // mark the parent for update and clear its conditionals if needed
-        {
             //printf ("Resetting parentCode = %ld\n", parentCode);
             taggedNodes.list_data[parentCode]     = 1L;
             hyFloat     const *localScalingFactor      = scalingAdjustments + parentCode*siteCount;
@@ -1616,10 +1664,11 @@ void            _TheTree::ComputeBranchCache    (
         _CalcNode    * currentTreeNode = (_CalcNode*) (isLeaf?  flatCLeaves (nodeCode):
                                                        flatTree    (notPassedRoot?nodeCode:parentCode));
         
-         hyFloat  const *  transitionMatrix = currentTreeNode->GetCompExp(catID)->theData;
-        
-
+        /*if (likeFuncEvalCallCount == 15098) {
+            fprintf (stderr, "%ld/%ld (%ld/%ld) => %s\n", nodeID, nodeCode, isLeaf, notPassedRoot, currentTreeNode->GetName()->get_str());
+        }*/
         
+        hyFloat  const *  transitionMatrix = currentTreeNode->GetCompExp(catID)->theData;
         hyFloat  *       childVector,*     lastUpdatedSite;
         
         if (!isLeaf) {
@@ -1630,8 +1679,7 @@ void            _TheTree::ComputeBranchCache    (
         long currentTCCIndex,currentTCCBit,parentTCCIIndex,parentTCCIBit;
         __ll_handle_tcc_init (tcc, isLeaf, siteCount, siteFrom, nodeCode, parentCode, parentTCCIBit, parentTCCIIndex, currentTCCBit, currentTCCIndex);
 
-        bool canScale = !notPassedRoot;
-
+ 
         if (alphabetDimension == 4L) {
             #ifdef _SLKP_USE_AVX_INTRINSICS
                 __m256d tmatrix_transpose [4] = {
@@ -1642,9 +1690,13 @@ void            _TheTree::ComputeBranchCache    (
                 };
             #endif
             for (long siteID = siteFrom; siteID < siteTo; siteID++, parentConditionals += 4L) {
+                bool canScale = !notPassedRoot;
                 hyFloat  const *tMatrix = transitionMatrix;
                 if (__lcache_loop_preface<4>(
                 isLeaf, lNodeFlags, siteID, siteOrdering, nodeCode, siteCount, siteFrom, parentConditionals, tMatrix, canScale, childVector, lastUpdatedSite, tcc, currentTCCBit, currentTCCIndex, parentTCCIBit, parentTCCIIndex, notPassedRoot, lNodeResolutions)) {
+                    /*if (likeFuncEvalCallCount == 15098 && siteID == 91) {
+                        fprintf (stderr, "__lcache_loop_preface (%ld) %g %g %g %g\n", nodeCode, parentConditionals[0], parentConditionals[1], parentConditionals[2], parentConditionals[3]);
+                    }*/
                     continue;
                 }
                 long     didScale =  0;
@@ -1659,6 +1711,9 @@ void            _TheTree::ComputeBranchCache    (
                     __ll_loop_handle_scaling<4L, false> (sum, parentConditionals, scalingAdjustments, didScale, nodeCode, siteCount, siteID, localScalerChange, theFilter->theFrequencies.get (siteOrdering.list_data[siteID]));
 
                 }
+                /*if (likeFuncEvalCallCount == 15098 && siteID == 91) {
+                    fprintf (stderr, "NODE = %ld, PARENT = %ld (%ld), P(G) = %lg, P(T) = %lg, scale = %ld\n", nodeCode, parentCode, canScale, parentConditionals[2], parentConditionals[3], didScale);
+                }*/
                 childVector += 4L;
                 __handle_site_corrections(didScale, siteID);
             }
@@ -1667,6 +1722,7 @@ void            _TheTree::ComputeBranchCache    (
                     __ll_handle_matrix_transpose<20>(transitionMatrix, tMatrixT);
                 #endif
                 for (long siteID = siteFrom; siteID < siteTo; siteID++, parentConditionals += 20L) {
+                    bool canScale = !notPassedRoot;
                     hyFloat  const *tMatrix = transitionMatrix;
                     if (__lcache_loop_preface<20>(
                     isLeaf, lNodeFlags, siteID, siteOrdering, nodeCode, siteCount, siteFrom, parentConditionals, tMatrix, canScale, childVector, lastUpdatedSite, tcc, currentTCCBit, currentTCCIndex, parentTCCIBit, parentTCCIIndex, notPassedRoot, lNodeResolutions)) {
@@ -1699,6 +1755,7 @@ void            _TheTree::ComputeBranchCache    (
                 __ll_handle_matrix_transpose<60L>(transitionMatrix, tMatrixT);
             #endif
             for (long siteID = siteFrom; siteID < siteTo; siteID++, parentConditionals += 60L) {
+                bool canScale = !notPassedRoot;
                 hyFloat  const *tMatrix = transitionMatrix;
                 if (__lcache_loop_preface<60>(
                 isLeaf, lNodeFlags, siteID, siteOrdering, nodeCode, siteCount, siteFrom, parentConditionals, tMatrix, canScale, childVector, lastUpdatedSite, tcc, currentTCCBit, currentTCCIndex, parentTCCIBit, parentTCCIIndex, notPassedRoot, lNodeResolutions)) {
@@ -1738,6 +1795,7 @@ void            _TheTree::ComputeBranchCache    (
                 __ll_handle_matrix_transpose<61L>(transitionMatrix, tMatrixT);
             #endif
             for (long siteID = siteFrom; siteID < siteTo; siteID++, parentConditionals += 61L) {
+                bool canScale = !notPassedRoot;
                 hyFloat  const *tMatrix = transitionMatrix;
                 if (__lcache_loop_preface<61>(
                 isLeaf, lNodeFlags, siteID, siteOrdering, nodeCode, siteCount, siteFrom, parentConditionals, tMatrix, canScale, childVector, lastUpdatedSite, tcc, currentTCCBit, currentTCCIndex, parentTCCIBit, parentTCCIIndex, notPassedRoot, lNodeResolutions)) {
@@ -1782,6 +1840,7 @@ void            _TheTree::ComputeBranchCache    (
                     __ll_product_sum_loop<61L> (tMatrix, childVector, parentConditionals, sum);
                 #endif
                 if (canScale) {
+                    bool canScale = !notPassedRoot;
                     #if defined _SLKP_USE_AVX_INTRINSICS
                         sum = _avx_sum_4(grandTotal) + s60;
                     #elif defined _SLKP_USE_SSE_INTRINSICS
@@ -1797,6 +1856,7 @@ void            _TheTree::ComputeBranchCache    (
                 __ll_handle_matrix_transpose<62L>(transitionMatrix, tMatrixT);
             #endif
             for (long siteID = siteFrom; siteID < siteTo; siteID++, parentConditionals += 62L) {
+                bool canScale = !notPassedRoot;
                 hyFloat  const *tMatrix = transitionMatrix;
                 if (__lcache_loop_preface<62>(
                 isLeaf, lNodeFlags, siteID, siteOrdering, nodeCode, siteCount, siteFrom, parentConditionals, tMatrix, canScale, childVector, lastUpdatedSite, tcc, currentTCCBit, currentTCCIndex, parentTCCIBit, parentTCCIIndex, notPassedRoot, lNodeResolutions)) {
@@ -1876,6 +1936,7 @@ void            _TheTree::ComputeBranchCache    (
                 __ll_handle_matrix_transpose<63L>(transitionMatrix, tMatrixT);
             #endif
             for (long siteID = siteFrom; siteID < siteTo; siteID++, parentConditionals += 63L) {
+                bool canScale = !notPassedRoot;
                 hyFloat  const *tMatrix = transitionMatrix;
                 if (__lcache_loop_preface<63>(
                 isLeaf, lNodeFlags, siteID, siteOrdering, nodeCode, siteCount, siteFrom, parentConditionals, tMatrix, canScale, childVector, lastUpdatedSite, tcc, currentTCCBit, currentTCCIndex, parentTCCIBit, parentTCCIIndex, notPassedRoot, lNodeResolutions)) {
@@ -1972,6 +2033,8 @@ void            _TheTree::ComputeBranchCache    (
         } else  {
               
               for (long siteID = siteFrom; siteID < siteTo; siteID++, parentConditionals += alphabetDimension) {
+                bool canScale = !notPassedRoot;
+
                 hyFloat  const *tMatrix = transitionMatrix;
                 if (__lcache_loop_preface_generic(
                 isLeaf, lNodeFlags, siteID, siteOrdering, nodeCode, siteCount, siteFrom, parentConditionals, tMatrix, canScale, childVector, lastUpdatedSite, tcc, currentTCCBit, currentTCCIndex, parentTCCIBit, parentTCCIIndex, notPassedRoot, lNodeResolutions, alphabetDimension)) {
@@ -2007,7 +2070,9 @@ void            _TheTree::ComputeBranchCache    (
     const unsigned long site_bound = alphabetDimension*siteTo;
     for (unsigned long ii = siteFrom * alphabetDimension; ii < site_bound; ii++) {
         state[ii] = rootConditionals[ii];
-        //printf ("Root conditional [%ld] = %g, node state [%ld] = %g\n", ii, state[ii], ii, cache[ii]);
+        /*if (likeFuncEvalCallCount == 15098 && ii / alphabetDimension == 91) {
+            printf ("Site %ld, Root conditional [%ld] = %g, node state [%ld] = %g\n", ii/alphabetDimension, ii, state[ii], ii, cache[ii]);
+        }*/
     }
     
     if (!siteCorrectionCounts && localScalerChange) {
diff --git a/tests/hbltests/Results/HIVSweden.out b/tests/hbltests/Results/HIVSweden.out
index 6270702a4..79b8035b5 100644
--- a/tests/hbltests/Results/HIVSweden.out
+++ b/tests/hbltests/Results/HIVSweden.out
@@ -2,27 +2,27 @@
 *** RUNNING SINGLE RATE MODEL ***
 #################################
 
->Done in 6 seconds 
+>Done in 5 seconds 
 
--1137.68842968525
+-1137.68842968708
 
 -----------------------------------
 
-dN/dS = 0.9051395111492303
+dN/dS = 0.9051395450601901
 
 
 *** RUNNING MODEL 1 (Neutral) ***
 ######################################
 
->Done in 4 seconds 
+>Done in 3 seconds 
 
--1114.64198020114
+-1114.64197979838
 
 ------------------------------------------------
 
-dN/dS = 0.5549192853457614 (sample variance = 0.2120268218125002)
+dN/dS = 0.5549192972885679 (sample variance = 0.2120268320358276)
 
-Rate[1]=  0.07854092 (weight=0.4830173)
+Rate[1]=  0.07854090 (weight=0.4830173)
 Rate[2]=  1.00000000 (weight=0.5169827)
 
 ------------------------------------------------
@@ -37,215 +37,215 @@ Import the following part into a data processing program
 for further analysis
 
 
-Rate/Site	Rate=0.07854092332840576	Rate=1
-1	0.0004154868030844993	0.9995845131969154
-2	0.9494802777599012	0.05051972224009882
-3	0.852714319138601	0.1472856808613991
-4	0.9772954463652723	0.02270455363472769
-5	0.7089276354636231	0.2910723645363769
-6	0.9706035417626925	0.02939645823730756
-7	0.8411768348369065	0.1588231651630935
-8	0.6448515276617085	0.3551484723382916
-9	0.009023917300298745	0.9909760826997013
-10	0.2515943370540133	0.7484056629459866
-11	0.936297677903	0.06370232209699996
-12	0.8039914370800333	0.1960085629199667
-13	0.9592993214137728	0.04070067858622721
-14	0.6831220160298361	0.3168779839701638
-15	0.9686855178831441	0.03131448211685599
-16	0.9686855178831441	0.03131448211685599
-17	0.9494802777599012	0.05051972224009882
-18	0.09627728498761287	0.9037227150123872
-19	0.8046581788776952	0.1953418211223048
-20	0.4254706661865902	0.5745293338134099
-21	0.001162124101353502	0.9988378758986465
-22	0.001180622306850325	0.9988193776931497
-23	0.9494802777599012	0.05051972224009882
-24	7.060704132933913e-05	0.9999293929586708
-25	0.8600358212119197	0.1399641787880803
-26	2.239769512514203e-05	0.9999776023048749
-27	0.7157867575924926	0.2842132424075075
-28	1.823975588109884e-08	0.9999999817602442
-29	0.9785833755496693	0.02141662445033066
-30	0.5823851349951137	0.4176148650048864
-31	0.002379375352658659	0.9976206246473412
-32	0.8411768348369065	0.1588231651630935
-33	0.936297677903	0.06370232209699996
-34	0.69798372452262	0.30201627547738
-35	0.7783934327640782	0.2216065672359217
-36	0.01421065602379084	0.9857893439762092
-37	0.06336725390958133	0.9366327460904187
-38	0.6912845829568349	0.3087154170431651
-39	0.001919864543239184	0.9980801354567608
-40	0.003960918088554947	0.996039081911445
-41	0.5294505098102081	0.4705494901897919
-42	0.8383799272270489	0.161620072772951
-43	0.5921773317964253	0.4078226682035747
-44	0.2081955268461745	0.7918044731538254
-45	0.4812556782883108	0.5187443217116893
-46	0.01854863190175804	0.981451368098242
-47	0.7638297673435148	0.2361702326564851
-48	0.140281306330095	0.8597186936699051
-49	0.07649047470366999	0.92350952529633
-50	0.5565561271083909	0.4434438728916091
-51	4.40144046907096e-06	0.9999955985595309
-52	0.9686855178831441	0.03131448211685599
-53	0.9337996664820075	0.06620033351799245
-54	0.05713139133742776	0.9428686086625723
-55	0.6733977936690815	0.3266022063309185
-56	0.9785833755496693	0.02141662445033066
-57	0.03639030672341932	0.9636096932765807
-58	0.9366003459987027	0.06339965400129738
-59	0.02363166727159621	0.9763683327284037
-60	0.7157867575924926	0.2842132424075075
-61	0.274064141330423	0.7259358586695771
-62	0.1720633735589035	0.8279366264410964
-63	0.0669606680318964	0.9330393319681036
-64	0.159026375951253	0.840973624048747
-65	0.1371565849170371	0.8628434150829629
-66	2.834615995562944e-08	0.9999999716538401
-67	0.5390463378705244	0.4609536621294756
-68	0.0006628202514785398	0.9993371797485214
-69	2.873458507277488e-05	0.9999712654149272
-70	0.3234742246568454	0.6765257753431545
-71	0.224929320914983	0.7750706790850169
-72	0.1619285389636656	0.8380714610363343
-73	0.1394769516923916	0.8605230483076084
-74	0.9494802777599012	0.05051972224009882
-75	0.02458208685945937	0.9754179131405406
-76	6.411330212760092e-05	0.9999358866978725
-77	0.9715443626189237	0.02845563738107615
-78	0.7862794565538135	0.2137205434461865
-79	0.7328751080906223	0.2671248919093778
-80	0.9706035417626925	0.02939645823730756
-81	0.8611716902269267	0.1388283097730733
-82	0.7668232725859022	0.2331767274140977
-83	5.67272509224703e-05	0.9999432727490775
-84	0.006741794211132465	0.9932582057888676
-85	0.9683654697043729	0.03163453029562719
-86	0.9686855178831441	0.03131448211685599
-87	1.53110154208541e-05	0.9999846889845792
-88	0.7668232725859022	0.2331767274140977
-89	0.9316122623318611	0.06838773766813891
-90	0.03861791813650284	0.9613820818634972
-91	0.5872085267662471	0.4127914732337529
+Rate/Site	Rate=0.07854089951899505	Rate=1
+1	0.0004154841702208897	0.9995845158297791
+2	0.9494800946553694	0.05051990534463052
+3	0.8527138076557068	0.1472861923442931
+4	0.9772953209605599	0.02270467903944003
+5	0.7089270358876215	0.2910729641123785
+6	0.9706034286518881	0.02939657134811178
+7	0.8411759549187992	0.1588240450812008
+8	0.644850909128687	0.3551490908713131
+9	0.009023903715709264	0.9909760962842907
+10	0.2515932185849544	0.7484067814150457
+11	0.9362971760831602	0.0637028239168398
+12	0.8039911310919982	0.1960088689080018
+13	0.9592990339188422	0.04070096608115781
+14	0.6831216537039989	0.3168783462960011
+15	0.9686853075822512	0.03131469241774875
+16	0.9686853075822512	0.03131469241774875
+17	0.9494800946553694	0.05051990534463052
+18	0.09627709694116433	0.9037229030588357
+19	0.8046571650681456	0.1953428349318545
+20	0.4254688557117581	0.5745311442882419
+21	0.001162113715308391	0.9988378862846916
+22	0.001180617896019509	0.9988193821039805
+23	0.9494800946553694	0.05051990534463052
+24	7.060663427640608e-05	0.9999293933657236
+25	0.8600353102504529	0.1399646897495471
+26	2.239754017472222e-05	0.9999776024598254
+27	0.7157862269098134	0.2842137730901866
+28	1.823968038733901e-08	0.9999999817603197
+29	0.9785832556654707	0.02141674433452926
+30	0.582384830553312	0.417615169446688
+31	0.002379358355164586	0.9976206416448354
+32	0.8411759549187992	0.1588240450812008
+33	0.9362971760831602	0.0637028239168398
+34	0.6979822928927298	0.3020177071072702
+35	0.778392401117516	0.2216075988824839
+36	0.01421053968210576	0.9857894603178943
+37	0.06336691422415654	0.9366330857758434
+38	0.6912831003120626	0.3087168996879373
+39	0.001919852660566921	0.998080147339433
+40	0.003960890204721402	0.9960391097952785
+41	0.5294498773319029	0.470550122668097
+42	0.8383792474897703	0.1616207525102296
+43	0.5921771900482152	0.4078228099517847
+44	0.208194488060306	0.791805511939694
+45	0.4812547391790605	0.5187452608209394
+46	0.01854854305069108	0.9814514569493089
+47	0.7638290038081602	0.2361709961918397
+48	0.1402805111076657	0.8597194888923344
+49	0.07649019370504373	0.9235098062949563
+50	0.5565553788637815	0.4434446211362185
+51	4.401412699981744e-06	0.9999955985873
+52	0.9686853075822512	0.03131469241774875
+53	0.9337995342113474	0.06620046578865253
+54	0.05713111117495385	0.9428688888250462
+55	0.6733962368122044	0.3266037631877956
+56	0.9785832556654707	0.02141674433452926
+57	0.03638999158091847	0.9636100084190816
+58	0.9366001361131814	0.06339986388681852
+59	0.02363153610583355	0.9763684638941664
+60	0.7157862269098134	0.2842137730901866
+61	0.2740628614189302	0.7259371385810699
+62	0.1720628516827533	0.8279371483172466
+63	0.06696027896282203	0.933039721037178
+64	0.1590255451412893	0.8409744548587106
+65	0.1371557240403739	0.8628442759596261
+66	2.83460436728884e-08	0.9999999716539563
+67	0.5390460493463879	0.4609539506536121
+68	0.0006628127730370201	0.999337187226963
+69	2.873441346176776e-05	0.9999712655865381
+70	0.3234733611365553	0.6765266388634448
+71	0.2249283193971887	0.7750716806028114
+72	0.161927560424059	0.8380724395759409
+73	0.1394764462745381	0.8605235537254619
+74	0.9494800946553694	0.05051990534463052
+75	0.02458202885075857	0.9754179711492413
+76	6.411295676480708e-05	0.9999358870432352
+77	0.9715441430178714	0.02845585698212867
+78	0.7862785460776877	0.2137214539223124
+79	0.7328737186780433	0.2671262813219566
+80	0.9706034286518881	0.02939657134811178
+81	0.8611706364190057	0.1388293635809944
+82	0.7668224471831266	0.2331775528168733
+83	5.672677285259407e-05	0.9999432732271474
+84	0.006741759929513022	0.993258240070487
+85	0.96836526712463	0.03163473287536989
+86	0.9686853075822512	0.03131469241774875
+87	1.531106965534881e-05	0.9999846889303446
+88	0.7668224471831266	0.2331775528168733
+89	0.9316117506366601	0.06838824936333981
+90	0.0386175672192274	0.9613824327807726
+91	0.5872082036930426	0.4127917963069574
 
 *** RUNNING MODEL 2 (Selection) ***
 ######################################
 
->Done in 11 seconds 
+>Done in 8 seconds 
 
--1106.44557162008
+-1106.44545333916
 
 ------------------------------------------------
 
-dN/dS = 1.130289110692444 (sample variance = 1.591177897940909)
+dN/dS = 1.129068005488161 (sample variance = 1.588067143267372)
 
-Rate[1]=  0.05803552 (weight=0.3738941)
-Rate[2]=  1.00000000 (weight=0.4437247)
-Rate[3]=  3.64547162 (weight=0.1823811)
+Rate[1]=  0.05854806 (weight=0.3747723)
+Rate[2]=  1.00000000 (weight=0.4427396)
+Rate[3]=  3.64070986 (weight=0.1824881)
 
 ------------------------------------------------
 
  Sites with dN/dS>1 (Posterior cutoff = 0.9)
 
-26 (0.9072435350390498)
-28 (0.9992039672038171)
-66 (0.9984916647170328)
-87 (0.9856690772721554)
+26 (0.9071235888236596)
+28 (0.999208561635355)
+66 (0.9984823131311625)
+87 (0.985639180825136)
 
 
 ------------------------------------------------
 
  Sites with dN/dS<=1 (Posterior cutoff = 0.9)
 
-1 (0.5051655529606731)
-2 (0.0004586047622872992)
-3 (0.005910658864265369)
-4 (9.646181948062124e-05)
-5 (0.03170339502414456)
-6 (0.0001508897972951693)
-7 (0.006938725368587095)
-8 (0.05670024471475955)
-9 (0.7225295530504859)
-10 (0.09638049255699262)
-11 (0.0009198237591565945)
-12 (0.01062846379830087)
-13 (0.000492759842270657)
-14 (0.03988301095601996)
-15 (0.0002222234934517179)
-16 (0.0002222234934517179)
-17 (0.0004586047622872992)
-18 (0.4327849002774984)
-19 (0.01160494123323598)
-20 (0.03363712511964109)
-21 (0.2470930495779124)
-22 (0.7984738153061429)
-23 (0.0004586047622872992)
-24 (0.5842524457215786)
-25 (0.005167728308291355)
-27 (0.02938020791912619)
-29 (7.937517626711708e-05)
-30 (0.09140259016602403)
-31 (0.5697739437706241)
-32 (0.006938725368587095)
-33 (0.0009198237591565945)
-34 (0.004510681205167178)
-35 (0.002049191785207676)
-36 (0.1050024130149768)
-37 (0.08580607126297962)
-38 (0.005393833634333682)
-39 (0.6439589837233456)
-40 (0.4558413160848798)
-41 (0.01183180050372747)
-42 (0.006852010985023359)
-43 (0.08361182186992157)
-44 (0.01179227732643408)
-45 (0.01748204559061378)
-46 (0.4290997677225385)
-47 (0.01836409760634284)
-48 (0.02526647483127544)
-49 (0.07094810710714877)
-50 (0.01061737774286497)
-51 (0.8862701608782172)
-52 (0.0002222234934517179)
-53 (0.0007670772046931027)
-54 (0.09812221578787755)
-55 (0.00633368543065168)
-56 (7.937517626711708e-05)
-57 (0.2048056654398318)
-58 (0.000726417634129623)
-59 (0.3311273587672101)
-60 (0.02938020791912619)
-61 (0.08189294413402423)
-62 (0.1934504467092352)
-63 (0.07706489107214247)
-64 (0.02140744778728731)
-65 (0.02736638268870059)
-67 (0.1261045684713601)
-68 (0.6055404804211916)
-69 (0.833597533314407)
-70 (0.05525842075956964)
-71 (0.1226377952204414)
-72 (0.02103774540341517)
-73 (0.2774126502695526)
-74 (0.0004586047622872992)
-75 (0.3069616959027847)
-76 (0.6757852477037506)
-77 (0.0002302410972905876)
-78 (0.01498418396816613)
-79 (0.003223456269673978)
-80 (0.0001508897972951693)
-81 (0.005171093611428542)
-82 (0.01775770256526435)
-83 (0.8143632663221626)
-84 (0.2774516014562403)
-85 (0.0001961741149948456)
-86 (0.0002222234934517179)
-88 (0.01775770256526435)
-89 (0.001172350305589507)
-90 (0.186172954156048)
-91 (0.08822774719946876)
+1 (0.5052373923971539)
+2 (0.0004593492021371499)
+3 (0.005917205452138728)
+4 (9.663607149571083e-05)
+5 (0.03172842746069128)
+6 (0.0001511795330813065)
+7 (0.006943635909993886)
+8 (0.05673891477006113)
+9 (0.7226692280451891)
+10 (0.09651028335914213)
+11 (0.0009207100403048468)
+12 (0.01064028334210353)
+13 (0.0004934182636844576)
+14 (0.03991745163515635)
+15 (0.0002225580972989534)
+16 (0.0002225580972989534)
+17 (0.0004593492021371499)
+18 (0.4330406113257519)
+19 (0.0116121595707726)
+20 (0.0336812685288552)
+21 (0.2473118054259542)
+22 (0.7984465921549734)
+23 (0.0004593492021371499)
+24 (0.5842280859182822)
+25 (0.005172980464989402)
+27 (0.02940476216359821)
+29 (7.952947195176931e-05)
+30 (0.09146670281887008)
+31 (0.5699030448214001)
+32 (0.006943635909993886)
+33 (0.0009207100403048468)
+34 (0.004515364181301996)
+35 (0.002050750658840959)
+36 (0.1052215418403245)
+37 (0.08602792425933428)
+38 (0.005398479347200226)
+39 (0.6440044083471241)
+40 (0.4558311806110331)
+41 (0.01185264798035508)
+42 (0.006857989522987185)
+43 (0.08367248274274289)
+44 (0.01183660486503293)
+45 (0.01750623812172373)
+46 (0.4294127609528174)
+47 (0.0183780006215666)
+48 (0.02534884977862542)
+49 (0.07112745644118525)
+50 (0.01063437293207849)
+51 (0.8861014582086777)
+52 (0.0002225580972989534)
+53 (0.0007683604393760527)
+54 (0.09835228649247904)
+55 (0.006339545564899368)
+56 (7.952947195176931e-05)
+57 (0.2050735349427733)
+58 (0.000727531168723411)
+59 (0.3314435529804894)
+60 (0.02940476216359821)
+61 (0.08200111674080049)
+62 (0.193666026549605)
+63 (0.0772513986346591)
+64 (0.02148421691961612)
+65 (0.02745286636506417)
+67 (0.1261823647518457)
+68 (0.605555060135807)
+69 (0.8333573458992974)
+70 (0.05534546347214257)
+71 (0.1227852244756305)
+72 (0.02110507208988902)
+73 (0.2776010249671293)
+74 (0.0004593492021371499)
+75 (0.3073445506286841)
+76 (0.6759845714674814)
+77 (0.0002305219770671781)
+78 (0.01499353952539858)
+79 (0.003225300542333516)
+80 (0.0001511795330813065)
+81 (0.005173815704594973)
+82 (0.01777066638393202)
+83 (0.8142466349701873)
+84 (0.2783391381741632)
+85 (0.0001964834081581603)
+86 (0.0002225580972989534)
+88 (0.01777066638393202)
+89 (0.001173394928906147)
+90 (0.1864176575163458)
+91 (0.08828981163248605)
 
 
 ------------------------------------------------
@@ -256,95 +256,95 @@ Import the following part into a data processing program
 for further analysis
 
 
-Rate/Site	Rate=0.05803551904096559	Rate=1	Rate=3.64547161773776
-1	2.039042780824998e-05	0.4948140566115184	0.5051655529606731
-2	0.8447543711358414	0.1547870241018713	0.0004586047622872992
-3	0.7161602740418139	0.2779290670939206	0.005910658864265369
-4	0.9018001736475444	0.09810336453297497	9.646181948062124e-05
-5	0.5817941052410429	0.3865024997348125	0.03170339502414456
-6	0.8869026108142981	0.1129464993884068	0.0001508897972951693
-7	0.7002133060545321	0.2928479685768806	0.006938725368587095
-8	0.5247416076674776	0.4185581476177629	0.05670024471475955
-9	0.001030497098954651	0.2764399498505593	0.7225295530504859
-10	0.09922289423382447	0.804396613209183	0.09638049255699262
-11	0.8185856564565643	0.1804945197842791	0.0009198237591565945
-12	0.6706899182118367	0.3186816179898624	0.01062846379830087
-13	0.8639581578245099	0.1355490823332193	0.000492759842270657
-14	0.5613226223127381	0.3987943667312419	0.03988301095601996
-15	0.8816214532222431	0.1181563232843053	0.0002222234934517179
-16	0.8816214532222431	0.1181563232843053	0.0002222234934517179
-17	0.8447543711358414	0.1547870241018713	0.0004586047622872992
-18	0.03258041875304914	0.5346346809694525	0.4327849002774984
-19	0.6634627593773434	0.3249322993894206	0.01160494123323598
-20	0.1578423652751871	0.8085205096051719	0.03363712511964109
-21	6.238970153028547e-05	0.7528445607205574	0.2470930495779124
-22	5.957427093987687e-05	0.2014666104229171	0.7984738153061429
-23	0.8447543711358414	0.1547870241018713	0.0004586047622872992
-24	1.722473813988775e-06	0.4157458318046073	0.5842524457215786
-25	0.7248831494482354	0.2699491222434732	0.005167728308291355
-26	1.79098661219942e-07	0.09275628586228903	0.9072435350390498
-27	0.5885447664274595	0.3820750256534142	0.02938020791912619
-28	7.739153851019017e-14	0.0007960327961054919	0.9992039672038171
-29	0.905003826260647	0.09491679856308587	7.937517626711708e-05
-30	0.470463461189876	0.4381339486440999	0.09140259016602403
-31	0.0001918709301161739	0.4300341852992596	0.5697739437706241
-32	0.7002133060545321	0.2928479685768806	0.006938725368587095
-33	0.8185856564565643	0.1804945197842791	0.0009198237591565945
-34	0.2930556524542092	0.7024336663406237	0.004510681205167178
-35	0.3526141222772732	0.6453366859375191	0.002049191785207676
-36	0.001268385541760359	0.8937292014432627	0.1050024130149768
-37	0.01154829414988456	0.9026456345871359	0.08580607126297962
-38	0.2911048433357678	0.7035013230298984	0.005393833634333682
-39	0.0001396861588575296	0.355901330117797	0.6439589837233456
-40	0.0002448465541435865	0.5439138373609765	0.4558413160848798
-41	0.2097066876421463	0.7784615118541262	0.01183180050372747
-42	0.6996255088940389	0.2935224801209377	0.006852010985023359
-43	0.4830787052206162	0.4333094729094621	0.08361182186992157
-44	0.02944754568248423	0.9587601769910816	0.01179227732643408
-45	0.1871451780503118	0.7953727763590743	0.01748204559061378
-46	0.003286918686078104	0.5676133135913832	0.4290997677225385
-47	0.6289631399679325	0.3526727624257246	0.01836409760634284
-48	0.02157260676327798	0.9531609184054467	0.02526647483127544
-49	0.01326287909287238	0.9157890137999789	0.07094810710714877
-50	0.2216072057730993	0.7677754164840356	0.01061737774286497
-51	2.348084455917174e-08	0.1137298156409383	0.8862701608782172
-52	0.8816214532222431	0.1181563232843053	0.0002222234934517179
-53	0.8203687814430313	0.1788641413522757	0.0007670772046931027
-54	0.01075878173450459	0.891119002477618	0.09812221578787755
-55	0.2804057512651671	0.7132605633041813	0.00633368543065168
-56	0.905003826260647	0.09491679856308587	7.937517626711708e-05
-57	0.006889612047585912	0.7883047225125823	0.2048056654398318
-58	0.8232664662614319	0.1760071161044386	0.000726417634129623
-59	0.004480976740308936	0.664391664492481	0.3311273587672101
-60	0.5885447664274595	0.3820750256534142	0.02938020791912619
-61	0.1064914550524415	0.8116156008135343	0.08189294413402423
-62	0.06777171429588846	0.7387778389948764	0.1934504467092352
-63	0.01214975847301941	0.9107853504548381	0.07706489107214247
-64	0.02337288930759019	0.9552196629051225	0.02140744778728731
-65	0.02103532247533653	0.951598294835963	0.02736638268870059
-66	2.487852398499066e-12	0.001508335280479298	0.9984916647170328
-67	0.4285637071264611	0.4453317244021789	0.1261045684713601
-68	2.969924066760578e-05	0.3944298203381409	0.6055404804211916
-69	4.031355426340145e-07	0.1664020635500504	0.833597533314407
-70	0.1273454140934302	0.8173961651470001	0.05525842075956964
-71	0.08804086160724707	0.7893213431723116	0.1226377952204414
-72	0.02443240135146625	0.9545298532451186	0.02103774540341517
-73	0.05234553025846337	0.6702418194719839	0.2774126502695526
-74	0.8447543711358414	0.1547870241018713	0.0004586047622872992
-75	0.004873253850134324	0.688165050247081	0.3069616959027847
-76	1.215507587476473e-06	0.324213536788662	0.6757852477037506
-77	0.8877821246335438	0.1119876342691656	0.0002302410972905876
-78	0.6474357779532641	0.3375800380785697	0.01498418396816613
-79	0.3289377286767906	0.6678388150535354	0.003223456269673978
-80	0.8869026108142981	0.1129464993884068	0.0001508897972951693
-81	0.7182778008443462	0.2765511055442251	0.005171093611428542
-82	0.6308413283365097	0.3514009690982259	0.01775770256526435
-83	9.028853568587064e-07	0.1856358307924804	0.8143632663221626
-84	9.578642651884824e-05	0.7224526121172409	0.2774516014562403
-85	0.88015102949012	0.1196527963948851	0.0001961741149948456
-86	0.8816214532222431	0.1181563232843053	0.0002222234934517179
-87	1.458685175204862e-08	0.01433090814099288	0.9856690772721554
-88	0.6308413283365097	0.3514009690982259	0.01775770256526435
-89	0.8109004627962602	0.1879271868981504	0.001172350305589507
-90	0.007317318166084692	0.8065097276778674	0.186172954156048
-91	0.4748123679645572	0.4369598848359741	0.08822774719946876
+Rate/Site	Rate=0.05854805862606752	Rate=1	Rate=3.640709855148721
+1	2.123821587431194e-05	0.4947413693869718	0.5052373923971539
+2	0.8454805150711112	0.1540601357267516	0.0004593492021371499
+3	0.7171925566791651	0.2768902378686959	0.005917205452138728
+4	0.9023145803022071	0.09758878362629728	9.663607149571083e-05
+5	0.5829266831612094	0.3853448893780992	0.03172842746069128
+6	0.8874735538558555	0.112375266611063	0.0001511795330813065
+7	0.701299413788691	0.2917569503013152	0.006943635909993886
+8	0.5258486059883327	0.4174124792416061	0.05673891477006113
+9	0.00105308734360597	0.276277684611205	0.7226692280451891
+10	0.1004481920361048	0.803041524604753	0.09651028335914213
+11	0.8194140603392319	0.1796652296204631	0.0009207100403048468
+12	0.6717772510719154	0.3175824655859811	0.01064028334210353
+13	0.8646165745664404	0.1348900071698753	0.0004934182636844576
+14	0.5624375756998519	0.3976449726649919	0.03991745163515635
+15	0.8822191864215847	0.1175582554811163	0.0002225580972989534
+16	0.8822191864215847	0.1175582554811163	0.0002225580972989534
+17	0.8454805150711112	0.1540601357267516	0.0004593492021371499
+18	0.032983308732931	0.533976079941317	0.4330406113257519
+19	0.6645919677056602	0.3237958727235673	0.0116121595707726
+20	0.1597162578236267	0.806602473647518	0.0336812685288552
+21	6.497561809983018e-05	0.752623218955946	0.2473118054259542
+22	6.147883799728457e-05	0.2014919290070293	0.7984465921549734
+23	0.8454805150711112	0.1540601357267516	0.0004593492021371499
+24	1.810178742010458e-06	0.4157701039029758	0.5842280859182822
+25	0.7259101630303384	0.2689168565046722	0.005172980464989402
+26	1.883724016460931e-07	0.09287622280393867	0.9071235888236596
+27	0.5896752507993037	0.3809199870370982	0.02940476216359821
+28	7.593295246048614e-14	0.0007914383645691474	0.999208561635355
+29	0.9055057132759483	0.09441475725209994	7.952947195176931e-05
+30	0.4715021983865639	0.437031098794566	0.09146670281887008
+31	0.0001979673775266633	0.4298989878010734	0.5699030448214001
+32	0.701299413788691	0.2917569503013152	0.006943635909993886
+33	0.8194140603392319	0.1796652296204631	0.0009207100403048468
+34	0.2960737907520586	0.6994108450666394	0.004515364181301996
+35	0.3559604560803941	0.641988793260765	0.002050750658840959
+36	0.00130942674258278	0.8934690314170927	0.1052215418403245
+37	0.01180978530471919	0.9021622904359465	0.08602792425933428
+38	0.2941041942604504	0.7004973263923493	0.005398479347200226
+39	0.000144144499790302	0.3558514471530856	0.6440044083471241
+40	0.0002529699651473707	0.5439158494238194	0.4558311806110331
+41	0.2120674596489594	0.7760798923706854	0.01185264798035508
+42	0.7007026627516525	0.2924393477253602	0.006857989522987185
+43	0.4841294396612101	0.432198077596047	0.08367248274274289
+44	0.0301228872990572	0.9580405078359099	0.01183660486503293
+45	0.1893248839788894	0.7931688778993868	0.01750623812172373
+46	0.003359794088751816	0.5672274449584306	0.4294127609528174
+47	0.63009891120623	0.3515230881722035	0.0183780006215666
+48	0.02206748041476434	0.9525836698066102	0.02534884977862542
+49	0.01356591438428317	0.9153066291745316	0.07112745644118525
+50	0.2240745276745924	0.7652910993933291	0.01063437293207849
+51	2.492645830161824e-08	0.1138985168648641	0.8861014582086777
+52	0.8822191864215847	0.1175582554811163	0.0002225580972989534
+53	0.8211660693180184	0.1780655702426055	0.0007683604393760527
+54	0.01100295212811388	0.8906447613794071	0.09835228649247904
+55	0.2833340811759865	0.7103263732591142	0.006339545564899368
+56	0.9055057132759483	0.09441475725209994	7.952947195176931e-05
+57	0.007045975624482142	0.7878804894327445	0.2050735349427733
+58	0.8240613622898397	0.1752111065414368	0.000727531168723411
+59	0.004581163540158042	0.6639752834793526	0.3314435529804894
+60	0.5896752507993037	0.3809199870370982	0.02940476216359821
+61	0.107806169745878	0.8101927135133216	0.08200111674080049
+62	0.06861334373069095	0.7377206297197041	0.193666026549605
+63	0.01242704838986371	0.9103215529754771	0.0772513986346591
+64	0.02390636440935269	0.9546094186710312	0.02148421691961612
+65	0.02151934533944069	0.9510277882954953	0.02745286636506417
+66	2.671166083534726e-12	0.001517686866166226	0.9984823131311625
+67	0.4295450384222473	0.4442725968259071	0.1261823647518457
+68	3.079998533959148e-05	0.3944141398788535	0.605555060135807
+69	4.243743084525438e-07	0.1666422297263941	0.8333573458992974
+70	0.1288881665397853	0.8157663699880722	0.05534546347214257
+71	0.08913564852994389	0.7880791269944257	0.1227852244756305
+72	0.02499602234463322	0.9538989055654777	0.02110507208988902
+73	0.05300541149090041	0.6693935635419703	0.2776010249671293
+74	0.8454805150711112	0.1540601357267516	0.0004593492021371499
+75	0.004981355134669857	0.687674094236646	0.3073445506286841
+76	1.275841252079727e-06	0.3240141526912666	0.6759845714674814
+77	0.8883601983297144	0.1114092796932183	0.0002305219770671781
+78	0.6485688357833093	0.3364376246912923	0.01499353952539858
+79	0.3321969903648257	0.6645777090928409	0.003225300542333516
+80	0.8874735538558555	0.112375266611063	0.0001511795330813065
+81	0.7193517468744175	0.2754744374209876	0.005173815704594973
+82	0.6319797015005985	0.3502496321154693	0.01777066638393202
+83	9.451297855017863e-07	0.1857524199000273	0.8142466349701873
+84	9.483309866743794e-05	0.7215660287271692	0.2783391381741632
+85	0.8807566499603399	0.119046866631502	0.0001964834081581603
+86	0.8822191864215847	0.1175582554811163	0.0002225580972989534
+87	1.512672292179054e-08	0.01436080404814105	0.985639180825136
+88	0.6319797015005985	0.3502496321154693	0.01777066638393202
+89	0.8117475469052963	0.1870790581657977	0.001173394928906147
+90	0.007483970311633406	0.8060983721720208	0.1864176575163458
+91	0.4758571278249646	0.4358530605425494	0.08828981163248605
diff --git a/tests/hbltests/Results/HIVSweden.out_MODEL_-1.nex b/tests/hbltests/Results/HIVSweden.out_MODEL_-1.nex
index c851cb79b..f5c923a96 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.9051395111492303;
-global kappa=0.4052422164401616;
+global c=0.9051395450601901;
+global kappa=0.4052422178470169;
 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.1932757871041837;
-givenTree.U68497.t=0.6442091652256334;
-givenTree.Node3.t=1.244091415592038;
-givenTree.U68501.t=1.072314201483783;
-givenTree.U68502.t=1.555369080763556;
-givenTree.U68503.t=1.573795705681594;
-givenTree.U68504.t=0.154501784041082;
-givenTree.U68506.t=0.5724439233580223;
-givenTree.Node15.t=0.4520764356621616;
-givenTree.U68505.t=0.2464768788725834;
-givenTree.Node14.t=0.2982771643141667;
-givenTree.U68507.t=0.4974362344804469;
-givenTree.Node13.t=0.2831885640218609;
-givenTree.U68508.t=1.61391993653481;
-givenTree.Node12.t=0.4115609436825465;
-givenTree.Node10.t=0.5209728609075241;
-givenTree.Node8.t=0.2739994031353166;
-givenTree.Node6.t=0.3910850665208743;
-givenTree.Node2.t=0.1922076567933305;
-givenTree.U68498.t=0.4102499915642663;
-givenTree.U68499.t=0.2311741941445719;
-givenTree.U68500.t=0.8119659493967246;
-givenTree.Node6_1.t=0.6788164927550273;
+givenTree.U68496.t=0.1932757790297635;
+givenTree.U68497.t=0.6442091360697517;
+givenTree.Node3.t=1.244091350196699;
+givenTree.U68501.t=1.072314145612458;
+givenTree.U68502.t=1.555368999549274;
+givenTree.U68503.t=1.573795636757062;
+givenTree.U68504.t=0.1545017778755618;
+givenTree.U68506.t=0.5724438957384254;
+givenTree.Node15.t=0.452076415208067;
+givenTree.U68505.t=0.2464768693858635;
+givenTree.Node14.t=0.2982771493845721;
+givenTree.U68507.t=0.4974360955659346;
+givenTree.Node13.t=0.2831885574728508;
+givenTree.U68508.t=1.613919861191977;
+givenTree.Node12.t=0.41156092487671;
+givenTree.Node10.t=0.5209728334125625;
+givenTree.Node8.t=0.2739993974710216;
+givenTree.Node6.t=0.3910850502827137;
+givenTree.Node2.t=0.1922077082389846;
+givenTree.U68498.t=0.4102499738247144;
+givenTree.U68499.t=0.2311741894617957;
+givenTree.U68500.t=0.8119659118420136;
+givenTree.Node6_1.t=0.678816471820221;
 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 e19d14c55..e9838a0e4 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.07854092332840576;
+global W=0.07854089951899505;
 W:<1;
-global P=0.4830173427363876;
+global P=0.483017317295036;
 P:<1;
-global kappa=0.3863939052192686;
+global kappa=0.3863915314909612;
 
 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.3014658105724962;
-givenTree.U68497.t=1.058937175876235;
-givenTree.Node3.t=2.081748321243699;
-givenTree.U68501.t=1.764331517307145;
-givenTree.U68502.t=2.673570235134937;
-givenTree.U68503.t=2.699658423145014;
-givenTree.U68504.t=0.243154063524979;
-givenTree.U68506.t=0.9038078207510426;
-givenTree.Node15.t=0.7153494550371319;
-givenTree.U68505.t=0.3956537768124779;
-givenTree.Node14.t=0.4993396138326727;
-givenTree.U68507.t=0.7737343409382251;
-givenTree.Node13.t=0.2771439596307402;
-givenTree.U68508.t=2.851718598450606;
-givenTree.Node12.t=0.747826310312856;
-givenTree.Node10.t=0.8232043891596889;
-givenTree.Node8.t=0.3093558783899246;
-givenTree.Node6.t=0.6402994000640408;
-givenTree.Node2.t=0.3042886856398082;
-givenTree.U68498.t=0.6387905520441879;
-givenTree.U68499.t=0.3594049973891654;
-givenTree.U68500.t=1.319977303229773;
-givenTree.Node6_1.t=1.138047087568419;
+givenTree.U68496.t=0.3014658876783538;
+givenTree.U68497.t=1.058937965054783;
+givenTree.Node3.t=2.08174984846752;
+givenTree.U68501.t=1.764333870286974;
+givenTree.U68502.t=2.673570567299918;
+givenTree.U68503.t=2.699658069558426;
+givenTree.U68504.t=0.2431541588198822;
+givenTree.U68506.t=0.9038084521967604;
+givenTree.Node15.t=0.7153498168088975;
+givenTree.U68505.t=0.3956539001874164;
+givenTree.Node14.t=0.4993402949067045;
+givenTree.U68507.t=0.7737346315034068;
+givenTree.Node13.t=0.277145717534158;
+givenTree.U68508.t=2.851718049074558;
+givenTree.Node12.t=0.7478267266145866;
+givenTree.Node10.t=0.8232048503168657;
+givenTree.Node8.t=0.309357412100375;
+givenTree.Node6.t=0.6402997546627597;
+givenTree.Node2.t=0.3042887273137687;
+givenTree.U68498.t=0.6387909451494354;
+givenTree.U68499.t=0.3594052835330087;
+givenTree.U68500.t=1.319979691408115;
+givenTree.Node6_1.t=1.13804766964687;
 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 710dbdae7..27c62b245 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.708705612029199;
+global P2=0.7081253811493516;
 P2:<1;
-global W_1=0.05803551904096559;
+global W_1=0.05854805862606752;
 W_1:<1;
-global W_2=3.64547161773776;
+global W_2=3.640709855148721;
 W_2:>1;
-global P1=0.3738941254448849;
+global P1=0.3747723001916182;
 P1:<1;
-global kappa=0.3595872901133234;
+global kappa=0.3597040016204291;
 
 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.1818762919281609;
-givenTree.U68497.t=0.6300551440548584;
-givenTree.Node3.t=1.242477534620252;
-givenTree.U68501.t=1.052718749766851;
-givenTree.U68502.t=1.745690425208843;
-givenTree.U68503.t=1.772661631401772;
-givenTree.U68504.t=0.1461819892700968;
-givenTree.U68506.t=0.5197050910198305;
-givenTree.Node15.t=0.4162511948548048;
-givenTree.U68505.t=0.2408046326551326;
-givenTree.Node14.t=0.275235108607599;
-givenTree.U68507.t=0.4766147197519663;
-givenTree.Node13.t=0.007124632417842378;
-givenTree.U68508.t=1.940263089445053;
-givenTree.Node12.t=0.5074844030410871;
-givenTree.Node10.t=0.5064503192192543;
-givenTree.Node8.t=0.1807141019606832;
-givenTree.Node6.t=0.3934552618450088;
-givenTree.Node2.t=0.1894290890321033;
-givenTree.U68498.t=0.333581631786244;
-givenTree.U68499.t=0.1874655318656705;
-givenTree.U68500.t=0.7992046138178572;
-givenTree.Node6_1.t=0.7479255472650597;
+givenTree.U68496.t=0.1819372645612003;
+givenTree.U68497.t=0.6305701510272947;
+givenTree.Node3.t=1.244260990073252;
+givenTree.U68501.t=1.05208595250357;
+givenTree.U68502.t=1.746445328113835;
+givenTree.U68503.t=1.774938469552289;
+givenTree.U68504.t=0.1463724860825859;
+givenTree.U68506.t=0.5200057981767465;
+givenTree.Node15.t=0.4166914000694648;
+givenTree.U68505.t=0.241055361539026;
+givenTree.Node14.t=0.2754520695076969;
+givenTree.U68507.t=0.4771824473110387;
+givenTree.Node13.t=0.006305605103123437;
+givenTree.U68508.t=1.944259154790765;
+givenTree.Node12.t=0.5081406979260248;
+givenTree.Node10.t=0.5074584441922946;
+givenTree.Node8.t=0.1808898503811754;
+givenTree.Node6.t=0.3939031724767187;
+givenTree.Node2.t=0.1896556393922262;
+givenTree.U68498.t=0.3337924593752155;
+givenTree.U68499.t=0.1875734352334583;
+givenTree.U68500.t=0.7998185516002427;
+givenTree.Node6_1.t=0.7492860788646263;
 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/libv3/MEME-partitioned.wbf b/tests/hbltests/libv3/MEME-partitioned.wbf
index 42793c14b..3aa80c1e3 100644
--- a/tests/hbltests/libv3/MEME-partitioned.wbf
+++ b/tests/hbltests/libv3/MEME-partitioned.wbf
@@ -20,7 +20,7 @@ fscanf ("data/CD2.nex.MEME.json","Raw",json);
 
 
 assert (check_value (
-    ((meme.json["fits"])["Unconstrained model"])["Log Likelihood"], -3466.56, 0.001), "Incorrect log-likelihood for the Global MG94xREV model");
+    ((meme.json["fits"])["Global MG94xREV"])["Log Likelihood"], -3466.57, 0.001), "Incorrect log-likelihood for the Global MG94xREV model");
     
 
 p_values   = (((meme.json["MLE"])["content"])["0"])[-1][6];
@@ -60,10 +60,9 @@ function confirm_site (site, p, dict, kind) {
     return FALSE;
 }
 
-utility.ForEachPair (p_values,"_index_", "_p_",
-                 "  
-                   confirm_site (_index_[0], _p_, test.expected_positives, 1))
-                ");
+for (_index_,_p_; in; p_values) {
+    confirm_site (_index_, _p_, test.expected_positives, 1);
+}
 
 assert (check_value (
     test.lrt_sum, 20.33, 0.05), "More than 5% difference in cumulative LRT for positively selected sites");