[med-svn] [Git][med-team/hyphy][master] 5 commits: routine-update: New upstream version
Nilesh Patra
gitlab at salsa.debian.org
Thu Oct 22 16:02:47 BST 2020
Nilesh Patra pushed to branch master at Debian Med / hyphy
Commits:
6ce7d810 by Nilesh Patra at 2020-10-22T20:14:26+05:30
routine-update: New upstream version
- - - - -
d8a6fecc by Nilesh Patra at 2020-10-22T20:14:27+05:30
New upstream version 2.5.20+dfsg
- - - - -
0912f873 by Nilesh Patra at 2020-10-22T20:14:39+05:30
Update upstream source from tag 'upstream/2.5.20+dfsg'
Update to upstream version '2.5.20+dfsg'
with Debian dir 1d345f1fbea12737020970729595867be9a76504
- - - - -
12c99283 by Nilesh Patra at 2020-10-22T20:17:57+05:30
Remove files from Files-Excluded that no longer shipped
- - - - -
f9cd333f by Nilesh Patra at 2020-10-22T20:31:28+05:30
Update changelog
- - - - -
27 changed files:
- CMakeLists.txt
- debian/changelog
- debian/copyright
- res/TemplateBatchFiles/SelectionAnalyses/FEL.bf
- res/TemplateBatchFiles/SelectionAnalyses/MEME.bf
- res/TemplateBatchFiles/SelectionAnalyses/SLAC.bf
- res/TemplateBatchFiles/SelectionAnalyses/contrast-fel.bf
- res/TemplateBatchFiles/SelectionAnalyses/modules/selection_lib.ibf
- res/TemplateBatchFiles/libv3/UtilityFunctions.bf
- res/TemplateBatchFiles/libv3/all-terms.bf
- res/TemplateBatchFiles/libv3/tasks/alignments.bf
- res/TemplateBatchFiles/libv3/tasks/ancestral.bf
- src/core/alignment.cpp
- src/core/batchlanruntime.cpp
- src/core/formula.cpp
- src/core/global_things.cpp
- src/core/include/matrix.h
- src/core/likefunc.cpp
- src/core/matrix.cpp
- src/core/operation.cpp
- src/core/tree.cpp
- src/core/tree_evaluator.cpp
- tests/hbltests/Results/HIVSweden.out
- tests/hbltests/Results/HIVSweden.out_MODEL_-1.nex
- tests/hbltests/Results/HIVSweden.out_MODEL_0.nex
- tests/hbltests/Results/HIVSweden.out_MODEL_1.nex
- tests/hbltests/libv3/MEME-partitioned.wbf
Changes:
=====================================
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)
=====================================
debian/changelog
=====================================
@@ -1,3 +1,10 @@
+hyphy (2.5.20+dfsg-1) unstable; urgency=medium
+
+ * Team upload.
+ * New upstream version
+
+ -- Nilesh Patra <npatra974 at gmail.com> Thu, 22 Oct 2020 20:14:59 +0530
+
hyphy (2.5.19+dfsg-1) unstable; urgency=medium
* New upstream version
=====================================
debian/copyright
=====================================
@@ -5,8 +5,6 @@ Source: https://github.com/veg/hyphy/releases
Files-Excluded: */SQLite*
.gitignore
*/.DS_Store
- */.xvpics
- */.gtkrc.swp
Files: *
Copyright: 1997-2015 Sergei L Kosakovsky Pond <spond at ucsd.edu>
=====================================
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,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);
}
+
=====================================
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 @@ private:
// matrix exp truncation precision
-
+ bool _validateCompressedStorage (void) const;
public:
@@ -518,9 +518,6 @@ public:
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;
=====================================
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;
}
//_____________________________________________________________________________________________
=====================================
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) {
=====================================
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) {
=====================================
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) {
=====================================
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
=====================================
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;
=====================================
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;
=====================================
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;
=====================================
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");
View it on GitLab: https://salsa.debian.org/med-team/hyphy/-/compare/b8c39452abc501bfc5ef6bee39675102ef402004...f9cd333fa9b7741889383e108de60ec3d4f98eb1
--
View it on GitLab: https://salsa.debian.org/med-team/hyphy/-/compare/b8c39452abc501bfc5ef6bee39675102ef402004...f9cd333fa9b7741889383e108de60ec3d4f98eb1
You're receiving this email because of your account on salsa.debian.org.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://alioth-lists.debian.net/pipermail/debian-med-commit/attachments/20201022/4de0684f/attachment-0001.html>
More information about the debian-med-commit
mailing list