Merge pull request #1235 from veg/develop

2.5.20

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)

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);

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,14 +279,17 @@ 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),
                                     Format(meme.report.row[3],7,3),
                                     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);
 }
+

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)

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];
            }

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",

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");
+            }
+    ');
 }
-*/
+

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";
     }

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;
         }
 

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];

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;
                     }
 

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);

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;
                 }
             }

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",

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) {

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;