[med-svn] [Git][med-team/hyphy][upstream] New upstream version 2.5.59+dfsg
Étienne Mollier (@emollier)
gitlab at salsa.debian.org
Sat Jan 13 16:07:13 GMT 2024
Étienne Mollier pushed to branch upstream at Debian Med / hyphy
Commits:
14ff91d5 by Étienne Mollier at 2024-01-13T16:59:21+01:00
New upstream version 2.5.59+dfsg
- - - - -
12 changed files:
- CMakeLists.txt
- res/TemplateBatchFiles/SelectionAnalyses/FADE.bf
- res/TemplateBatchFiles/SelectionAnalyses/MEME.bf
- res/TemplateBatchFiles/SelectionAnalyses/error-filter.bf
- src/core/formula.cpp
- src/core/global_things.cpp
- src/core/include/formula.h
- src/core/include/matrix.h
- src/core/include/simplelist.h
- src/core/likefunc.cpp
- src/core/matrix.cpp
- src/core/tree.cpp
Changes:
=====================================
CMakeLists.txt
=====================================
@@ -169,7 +169,11 @@ set_source_files_properties(${SRC_CORE} ${SRC_NEW} {SRC_UTILS} PROPERTIES COMPIL
set(DEFAULT_WARNING_FLAGS " -w -Weffc++ -Wextra -Wall ")
set(DEFAULT_DEBUG_WARNING_FLAGS "-Wall -Wno-int-to-pointer-cast -Wno-conversion-null -Wno-sign-compare")
+set(ADDITIONAL_FLAGS "")
+if (DEBUGFLAGS)
+ set(ADDITIONAL_FLAGS "-g")
+endif()
if(CMAKE_COMPILER_IS_GNUCC OR CMAKE_COMPILER_IS_GNUCXX)
@@ -289,7 +293,7 @@ if (CMAKE_CXX_COMPILER_ID MATCHES "Clang")
endif (CMAKE_CXX_COMPILER_ID MATCHES "Clang")
#set(DEFAULT_COMPILE_FLAGS "${DEFAULT_COMPILE_FLAGS} -fopt-info -fopt-info-vec-missed")
-set(DEFAULT_COMPILE_FLAGS "${DEFAULT_COMPILE_FLAGS}")
+set(DEFAULT_COMPILE_FLAGS "${DEFAULT_COMPILE_FLAGS} ${ADDITIONAL_FLAGS}")
if(NOT DEFINED DEFAULT_COMPILE_FLAGS)
set(DEFAULT_COMPILE_FLAGS "")
@@ -299,7 +303,7 @@ if(NOT DEFINED DEFAULT_LINK_FLAGS)
set(DEFAULT_LINK_FLAGS "")
endif(NOT DEFINED DEFAULT_LINK_FLAGS)
-set(DEFAULT_LINK_FLAGS "${DEFAULT_LINK_FLAGS} ")
+set(DEFAULT_LINK_FLAGS "${DEFAULT_LINK_FLAGS} ${ADDITIONAL_FLAGS}")
if(NOT DEFINED DEFAULT_WARNING_FLAGS)
set(DEFAULT_WARNING_FLAGS "")
@@ -309,7 +313,8 @@ if(CMAKE_SYSTEM_NAME STREQUAL "Emscripten")
set(DEFAULT_COMPILE_FLAGS "${DEFAULT_COMPILE_FLAGS} -fwasm-exceptions ")
endif()
-MESSAGE ("Set default compiler flags to ${DEFAULT_COMPILE_FLAGS}")
+MESSAGE ("Set compiler flags to ${DEFAULT_COMPILE_FLAGS}")
+MESSAGE ("Set linker flags to ${DEFAULT_LINK_FLAGS}")
#-------------------------------------------------------------------------------
# OpenMP support
=====================================
res/TemplateBatchFiles/SelectionAnalyses/FADE.bf
=====================================
@@ -544,7 +544,7 @@ for (fade.residue = 0; fade.residue < 20; fade.residue += 1) {
LikelihoodFunction fade.lf = (fade.lf.components);
estimators.ApplyExistingEstimates ("fade.lf", fade.model_id_to_object, fade.baseline_fit, None);
-
+
fade.conditionals.raw = fade.ComputeOnGrid ("fade.lf",
fade.grid.MatrixToDict (fade.cache[terms.fade.cache.grid]),
"fade.pass2.evaluator",
=====================================
res/TemplateBatchFiles/SelectionAnalyses/MEME.bf
=====================================
@@ -659,7 +659,15 @@ for (meme.partition_index = 0; meme.partition_index < meme.partition_count; meme
SetParameter (DEFER_CONSTRAINT_APPLICATION, 1, 0);
+ utility.ToggleEnvVariable ("LOOP_TREE_ITERATOR_PREORDER", TRUE);
+
+ is_root = TRUE;
+
for (_node_; in; meme.site_tree_fel) {
+ if (is_root) {
+ is_root = FALSE;
+ continue;
+ }
_bl_ = ((meme.final_partitioned_mg_results[terms.branch_length])[meme.partition_index])[_node_];
_node_class_ = (meme.selected_branches[meme.partition_index])[_node_];
if (_node_class_ != terms.tree_attributes.test) {
@@ -669,10 +677,10 @@ for (meme.partition_index = 0; meme.partition_index < meme.partition_count; meme
meme.apply_site_constraints ("meme.site_tree_bsrel",_node_,_bl_,meme.scaler_mapping ['FG']);
meme.apply_site_constraints ("meme.site_tree_fel",_node_,_bl_,meme.scaler_mapping ['FEL-FG']);
}
-
-
}
+ utility.ToggleEnvVariable ("LOOP_TREE_ITERATOR_PREORDER", None);
+
SetParameter (DEFER_CONSTRAINT_APPLICATION, 0, 0);
=====================================
res/TemplateBatchFiles/SelectionAnalyses/error-filter.bf
=====================================
@@ -94,10 +94,14 @@ if (efilter.error_sink_proportion == 0) {
}
efilter.error_sink_prior_odds_ratio = Min (1e25,efilter.error_sink_proportion / efilter.fast_omega_proportion);
+efilter.verbose_logging = utility.GetEnvVariable ("VERBOSE_LOGGING");
+
efilter.site_offset = 0;
utility.ToggleEnvVariable ("LOOP_TREE_ITERATOR_PREORDER", TRUE);
+
+
for (p = 0; p < efilter.input[terms.json.partition_count]; p+=1) {
efilter.sequences = {};
efilter.masked_sites = {};
@@ -175,6 +179,19 @@ for (p = 0; p < efilter.input[terms.json.partition_count]; p+=1) {
efilter.site_BF2 = 1e25;
}
+ if (efilter.verbose_logging ) {
+ if (efilter.site_BF >= efilter.threshold || efilter.site_BF2 >= efilter.ratio_threshold) {
+ if (efilter.site_BF >= efilter.threshold && efilter.site_BF2 >= efilter.ratio_threshold) {
+ fprintf (stdout, "FILTERED\n");
+ } else {
+ fprintf (stdout, "KEPT\n");
+ }
+ fprintf (stdout, ">" + (site+1) + "/" + node, " " + efilter.site_BF + ":" +efilter.site_BF2 + "\n");
+
+ }
+ }
+
+
if (efilter.site_BF >= efilter.threshold && efilter.site_BF2 >= efilter.ratio_threshold) {
if (efilter.masked_sites/node) { // terminal node
if (efilter.masked_already [ntm] != 1) {
@@ -189,12 +206,15 @@ for (p = 0; p < efilter.input[terms.json.partition_count]; p+=1) {
if (efilter.masked_already [ntm] != 1) {
efilter.masked_sites [ntm] + (site + efilter.site_offset);
}
- }
- //console.log ("Masking everything " + site + " " + node + " " + Abs (efilter.leaf_descendants) / Abs (efilter.sequences));
+ }
+ if (efilter.verbose_logging) {}
+ console.log ("Masking everything " + site + " " + node + " " + Abs (efilter.leaf_descendants) / Abs (efilter.sequences));
+ }
break;
}
for (ntm, ignore; in; efilter.leaf_descendants[node]) {
efilter.write_out[ntm] = "---";
+ //console.log ("Masking child " + ntm + " of parent " + node + " at site " + (1+site));
if (efilter.masked_already [ntm] != 1) {
efilter.masked_sites [ntm] + (site + efilter.site_offset);
}
=====================================
src/core/formula.cpp
=====================================
@@ -2309,15 +2309,26 @@ bool _Formula::ConvertToSimple (_AVLList& variable_index) {
delete [] simpleExpressionStatus;
simpleExpressionStatus = nil;
}
+
+ theStack.Reset();
+ long available_slots = MIN (32, sizeof (long) * MEMORYSTEP / sizeof (hyFloat));
+ long used_float_slots = 0;
+ hyFloat * store_constant = (hyFloat*)theStack.theStack._getStatic();
+
if (!theFormula.empty()) {
simpleExpressionStatus = new long [theFormula.countitems()];
for (unsigned long i=0UL; i<theFormula.countitems(); i++) {
_Operation* this_op = ItemAt (i);
- simpleExpressionStatus[i] = -4L;
+ simpleExpressionStatus[i] = -64L;
if (this_op->theNumber) {
- simpleExpressionStatus[i] = -1L;
+ if (used_float_slots < available_slots) {
+ store_constant [used_float_slots++] = this_op->theNumber->Value();
+ simpleExpressionStatus[i] = -1L - used_float_slots;
+ } else {
+ simpleExpressionStatus[i] = -1L;
+ }
} else if (this_op->theData >= 0) {
this_op->theData = variable_index.FindLong (this_op->theData);
simpleExpressionStatus[i] = this_op->theData;
@@ -2331,8 +2342,10 @@ bool _Formula::ConvertToSimple (_AVLList& variable_index) {
this_op->numberOfTerms = -3;
}
}
- if (this_op->opCode == HY_OP_CODE_RANDOM || this_op->opCode == HY_OP_CODE_TIME)
+
+ if (this_op->opCode == HY_OP_CODE_RANDOM || this_op->opCode == HY_OP_CODE_TIME) {
has_volatiles = true;
+ }
if (this_op->numberOfTerms == 2) {
switch (this_op->opCode) {
@@ -2393,95 +2406,98 @@ hyFloat _Formula::ComputeSimple (_SimpleFormulaDatum* stack, _SimpleFormulaDatum
if (theFormula.nonempty()) {
long stackTop = 0;
unsigned long upper_bound = NumberOperations();
-
+
+ const hyFloat * constants = (hyFloat*)theStack.theStack._getStatic();
+
for (unsigned long i=0UL; i<upper_bound; i++) {
- if (simpleExpressionStatus[i] == -1L) {
- stack[stackTop++].value = ItemAt (i)->theNumber->Value();
- /*if (thisOp->theNumber) {
- stack[stackTop++].value = thisOp->theNumber->Value();
- continue;*/
+ if (simpleExpressionStatus[i] >= 0L) {
+ stack[stackTop++] = varValues[simpleExpressionStatus[i]];
} else {
- if (simpleExpressionStatus[i] >= 0) {
- stack[stackTop++] = varValues[simpleExpressionStatus[i]];
+ if (simpleExpressionStatus[i] <= -2L && simpleExpressionStatus[i] >= -32L) {
+ stack[stackTop++].value = constants[-simpleExpressionStatus[i] - 2L];
} else {
- if (simpleExpressionStatus[i] <= -10000L) {
- stackTop--;
- switch (-simpleExpressionStatus[i] - 10000L) {
- case HY_OP_CODE_ADD:
- stack[stackTop-1].value = stack[stackTop-1].value + stack[stackTop].value;
- break;
- case HY_OP_CODE_SUB:
- stack[stackTop-1].value = stack[stackTop-1].value - stack[stackTop].value;
- break;
- case HY_OP_CODE_MUL:
- stack[stackTop-1].value = stack[stackTop-1].value * stack[stackTop].value;
- break;
- case HY_OP_CODE_DIV:
- stack[stackTop-1].value = stack[stackTop-1].value / stack[stackTop].value;
- break;
- case HY_OP_CODE_POWER: {
- //stack[stackTop-1].value = pow (stack[stackTop-1].value, stack[stackTop].value);
-
- if (stack[stackTop-1].value==0.0) {
- if (stack[stackTop].value > 0.0) {
- return 0.0;
- } else {
- return 1.0;
- }
- }
- stack[stackTop-1].value = pow (stack[stackTop-1].value, stack[stackTop].value);
-
- break;
- }
- default:
- HandleApplicationError ("Internal error in _Formula::ComputeSimple - unsupported shortcut operation.)", true);
- return 0.0;
- }
+ if (simpleExpressionStatus[i] == -1L) {
+ stack[stackTop++].value = ((_Operation**)theFormula.list_data)[i]->theNumber->Value();
} else {
- _Operation const* thisOp = ItemAt (i);
- stackTop--;
- if (thisOp->numberOfTerms==2) {
- hyFloat (*theFunc) (hyFloat, hyFloat);
- theFunc = (hyFloat(*)(hyFloat,hyFloat))thisOp->opCode;
- if (stackTop<1L) {
- HandleApplicationError ("Internal error in _Formula::ComputeSimple - stack underflow.)", true);
- return 0.0;
+ if (simpleExpressionStatus[i] <= -10000L) {
+ stackTop--;
+ switch (-simpleExpressionStatus[i] - 10000L) {
+ case HY_OP_CODE_ADD:
+ stack[stackTop-1].value = stack[stackTop-1].value + stack[stackTop].value;
+ break;
+ case HY_OP_CODE_SUB:
+ stack[stackTop-1].value = stack[stackTop-1].value - stack[stackTop].value;
+ break;
+ case HY_OP_CODE_MUL:
+ stack[stackTop-1].value = stack[stackTop-1].value * stack[stackTop].value;
+ break;
+ case HY_OP_CODE_DIV:
+ stack[stackTop-1].value = stack[stackTop-1].value / stack[stackTop].value;
+ break;
+ case HY_OP_CODE_POWER: {
+ //stack[stackTop-1].value = pow (stack[stackTop-1].value, stack[stackTop].value);
+
+ if (stack[stackTop-1].value==0.0) {
+ if (stack[stackTop].value > 0.0) {
+ return 0.0;
+ } else {
+ return 1.0;
+ }
+ }
+ stack[stackTop-1].value = pow (stack[stackTop-1].value, stack[stackTop].value);
+
+ break;
+ }
+ default:
+ HandleApplicationError ("Internal error in _Formula::ComputeSimple - unsupported shortcut operation.)", true);
+ return 0.0;
}
- stack[stackTop-1].value = (*theFunc)(stack[stackTop-1].value,stack[stackTop].value);
} else {
- switch (thisOp->numberOfTerms) {
- case -2 : {
- hyFloat (*theFunc) (hyPointer,hyFloat);
- theFunc = (hyFloat(*)(hyPointer,hyFloat))thisOp->opCode;
+ _Operation const* thisOp = ItemAt (i);
+ stackTop--;
+ if (thisOp->numberOfTerms==2) {
+ hyFloat (*theFunc) (hyFloat, hyFloat);
+ theFunc = (hyFloat(*)(hyFloat,hyFloat))thisOp->opCode;
if (stackTop<1L) {
HandleApplicationError ("Internal error in _Formula::ComputeSimple - stack underflow.)", true);
return 0.0;
}
- stack[stackTop-1].value = (*theFunc)(stack[stackTop-1].reference,stack[stackTop].value);
- break;
- }
- case -3 : {
- void (*theFunc) (hyPointer,hyFloat,hyFloat);
- theFunc = (void(*)(hyPointer,hyFloat,hyFloat))thisOp->opCode;
- if (stackTop != 2 || i != theFormula.lLength - 1) {
- HandleApplicationError ("Internal error in _Formula::ComputeSimple - stack underflow or MCoord command is not the last one.)", true);
-
- return 0.0;
- }
- //stackTop = 0;
- // value, reference, index
- (*theFunc)(stack[1].reference,stack[2].value, stack[0].value);
- break;
- }
- default: {
- hyFloat (*theFunc) (hyFloat);
- theFunc = (hyFloat(*)(hyFloat))thisOp->opCode;
- stack[stackTop].value = (*theFunc)(stack[stackTop].value);
- ++stackTop;
+ stack[stackTop-1].value = (*theFunc)(stack[stackTop-1].value,stack[stackTop].value);
+ } else {
+ switch (thisOp->numberOfTerms) {
+ case -2 : {
+ hyFloat (*theFunc) (hyPointer,hyFloat);
+ theFunc = (hyFloat(*)(hyPointer,hyFloat))thisOp->opCode;
+ if (stackTop<1L) {
+ HandleApplicationError ("Internal error in _Formula::ComputeSimple - stack underflow.)", true);
+ return 0.0;
+ }
+ stack[stackTop-1].value = (*theFunc)(stack[stackTop-1].reference,stack[stackTop].value);
+ break;
+ }
+ case -3 : {
+ void (*theFunc) (hyPointer,hyFloat,hyFloat);
+ theFunc = (void(*)(hyPointer,hyFloat,hyFloat))thisOp->opCode;
+ if (stackTop != 2 || i != theFormula.lLength - 1) {
+ HandleApplicationError ("Internal error in _Formula::ComputeSimple - stack underflow or MCoord command is not the last one.)", true);
+
+ return 0.0;
+ }
+ //stackTop = 0;
+ // value, reference, index
+ (*theFunc)(stack[1].reference,stack[2].value, stack[0].value);
+ break;
+ }
+ default: {
+ hyFloat (*theFunc) (hyFloat);
+ theFunc = (hyFloat(*)(hyFloat))thisOp->opCode;
+ stack[stackTop].value = (*theFunc)(stack[stackTop].value);
+ ++stackTop;
+ }
+ }
+
}
- }
-
}
}
}
=====================================
src/core/global_things.cpp
=====================================
@@ -122,7 +122,7 @@ namespace hy_global {
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"),
kErrorNumerical ("To treat numerical errors as warnings, please specify \"ENV=TOLERATE_NUMERICAL_ERRORS=1;\" as the command line argument. This often resolves the issue, which is indicative of numerical instability."),
- kHyPhyVersion = _String ("2.5.58"),
+ kHyPhyVersion = _String ("2.5.59"),
kNoneToken = "None",
kNullToken = "null",
=====================================
src/core/include/formula.h
=====================================
@@ -80,7 +80,7 @@ protected:
_List* resultCache;
_Stack theStack;
_List theFormula;
- long* simpleExpressionStatus;
+ long* simpleExpressionStatus;
/**
SLKP: 20200924
Added this shorthand to improve memory locality and speed-up SimpleCompute performance
=====================================
src/core/include/matrix.h
=====================================
@@ -75,6 +75,7 @@ struct _CompiledMatrixData {
hyFloat * formulaValues;
long * formulaRefs;
+ long stackDepth;
bool has_volatile_entries;
_SimpleList varIndex,
=====================================
src/core/include/simplelist.h
=====================================
@@ -684,6 +684,10 @@ class _SimpleList:public BaseObj {
long* quickArrayAccess(void) {
return (long*)list_data;
}
+
+ long* _getStatic () {
+ return static_data;
+ }
};
//TODO:Why is this a global function? If it needs to be, should be in helpers.cpp
=====================================
src/core/likefunc.cpp
=====================================
@@ -4589,10 +4589,13 @@ _Matrix* _LikelihoodFunction::Optimize (_AssociativeList const * options)
shrink_factor = MAX(GOLDEN_RATIO,MIN(stdFactor,oldAverage/averageChange));
long steps = logLHistory.get_used();
+ //printf ("DIFFS\n");
for (long k = 1; k <= MIN(5, steps-1); k++) {
diffs[k-1] = logLHistory.theData[steps-k] - logLHistory.theData[steps-k-1];
- //printf ("\n==> DIFFS %ld : %g\n", k, diffs[k-1]);
+ //printf ("==> DIFFS %ld : %g\n", k, diffs[k-1]);
}
+ //printf ("\n");
+
if (steps > 2 && diffs[0] >= diffs[1]) {
convergenceMode = 1;
}
@@ -4647,7 +4650,8 @@ _Matrix* _LikelihoodFunction::Optimize (_AssociativeList const * options)
}
}
- if (convergenceMode >= 2) {
+ //printf ("NO CHANGE %d\n", noChange.countitems());
+ if (convergenceMode >= 3) {
noChange.Clear();
noChange << -1;
}
@@ -4772,6 +4776,9 @@ _Matrix* _LikelihoodFunction::Optimize (_AssociativeList const * options)
}
large_change.Clear();
+
+ //ObjectToConsole (&_variables_that_dont_change);
+ //NLToConsole();
LoggerAddCoordinatewisePhase (shrink_factor, convergenceMode);
@@ -4810,12 +4817,20 @@ _Matrix* _LikelihoodFunction::Optimize (_AssociativeList const * options)
}
}
+ if (_variables_that_dont_change.get (current_index) > 1) {
+ if (convergenceMode <= 2 && inCount < termFactor -1) {
+ averageChange2+=bestVal;
+ continue;
+ }
+ }
+
_Vector *vH = nil;
hyFloat precisionStep = 0.,
brackStep;
if (use_adaptive_step) {
+
vH = (_Vector*)(*stepHistory)(current_index);
//hyFloat suggestedPrecision = currentPrecision*(1.+198./(1.+exp(sqrt(loopCounter))));
@@ -4934,9 +4949,17 @@ _Matrix* _LikelihoodFunction::Optimize (_AssociativeList const * options)
if (cj != 0.) {
averageChange += fabs (ch/cj);
}
- if ((ch < precisionStep*0.01 /*|| lastLogL - maxSoFar < 1e-6*/) && inCount == 0) {
- nc2 << current_index;
+ //printf ("VAR %d change %g (%g, %d)\n", current_index, ch, precisionStep*0.01, inCount);
+
+ if (ch < 1e-6) {
_variables_that_dont_change[current_index] ++;
+ } else {
+ _variables_that_dont_change[current_index] = 0;
+ }
+
+ if ((ch < precisionStep*0.01 /*|| lastLogL - maxSoFar < 1e-6*/) && inCount < termFactor - 1) {
+ nc2 << current_index;
+
}
} else {
averageChange += ch;
@@ -4999,7 +5022,15 @@ _Matrix* _LikelihoodFunction::Optimize (_AssociativeList const * options)
averageChange/=(hyFloat)(indexInd.lLength-nc2.lLength+1);
// mean of change in parameter values during the last step
+ //BufferToConsole("\nNO CHANGE BEFORE\n");
+ //ObjectToConsole(&noChange);
+ //NLToConsole ();
+ //ObjectToConsole(&nc2);
+ //NLToConsole ();
+
+
if (glVars.lLength == indexInd.lLength) {
+ // only global variables in the LF
noChange.Clear();
noChange.Duplicate (&nc2);
} else {
@@ -5015,20 +5046,24 @@ _Matrix* _LikelihoodFunction::Optimize (_AssociativeList const * options)
});
if (at_boundary.nonempty()) {
- _SimpleList nc3;
- nc3.Union (nc2,at_boundary);
- noChange.Subtract (nc3,glVars);
- //ObjectToConsole(&at_boundary);
- //NLToConsole();
+ //_SimpleList nc3;
+ //nc3.Union (nc2,at_boundary);
+ //noChange.Subtract (nc3,glVars);
+ noChange.Union (nc2, at_boundary);
+ // variables on the boundary and those that did not change will be skipped in the next iteration.
} else {
- noChange.Subtract (nc2,glVars);
+ //noChange.Subtract (nc2,glVars);
+ noChange = nc2;
}
+ // always interate over globals
}
if (noChange.empty()) {
noChange << -1;
}
-
+
+ //BufferToConsole("\nNO CHANGE AFTER\n");
+ //ObjectToConsole(&noChange);
forward = true;
=====================================
src/core/matrix.cpp
=====================================
@@ -2479,6 +2479,25 @@ bool _Matrix::ProcessFormulas (long& stackLength, _AVLList& varList, _S
_Formula ** theFormulas = (_Formula**)theData;
bool isGood = true;
+
+ auto process_formula = [&] (_Formula * this_formula) -> bool {
+ if (runAll || this_formula->AmISimple(stackLength,varList)) {
+ _String * flaString = (_String*)this_formula->toStr(kFormulaStringConversionNormal, nil,true);
+ long fref = flaStrings.Insert(flaString,newFormulas.lLength);
+ if (fref < 0) {
+ references << flaStrings.GetXtra (-fref-1);
+ DeleteObject (flaString);
+ } else {
+ newFormulas << (long)this_formula;
+ references << fref;
+ }
+
+ } else {
+ isGood = false;
+ return false;
+ }
+ return true;
+ };
if (theIndex) {
for (long i = 0L; i<lDim; i++) {
@@ -2488,20 +2507,7 @@ bool _Matrix::ProcessFormulas (long& stackLength, _AVLList& varList, _S
references << -1;
continue;
}
- thisFormula = theFormulas[i];
- if (runAll || thisFormula->AmISimple(stackLength,varList)) {
- _String * flaString = (_String*)thisFormula->toStr(kFormulaStringConversionNormal, nil,true);
- long fref = flaStrings.Insert(flaString,newFormulas.lLength);
- if (fref < 0) {
- references << flaStrings.GetXtra (-fref-1);
- DeleteObject (flaString);
- } else {
- newFormulas << (long)thisFormula;
- references << fref;
- }
-
- } else {
- isGood = false;
+ if (! process_formula (theFormulas[i])) {
break;
}
} else {
@@ -2510,28 +2516,18 @@ bool _Matrix::ProcessFormulas (long& stackLength, _AVLList& varList, _S
}
} else {
for (long i = 0L; i<lDim; i++) {
- if ((theFormulas[i]!=(_Formula*)ZEROPOINTER)&&(!theFormulas[i]->IsEmpty())) {
- thisFormula = theFormulas[i];
-
+
+ thisFormula = theFormulas[i];
+
+ if ((thisFormula!=(_Formula*)ZEROPOINTER)&&(!thisFormula->IsEmpty())) {
if (stencil && CheckEqual(stencil->theData[i],0.0)) {
references << -1;
continue;
}
-
- if (runAll || thisFormula->AmISimple(stackLength,varList)) {
- _String * flaString = (_String*)thisFormula->toStr(kFormulaStringConversionNormal, nil,true);
- long fref = flaStrings.Insert(flaString,newFormulas.lLength);
- if (fref < 0) {
- references << flaStrings.GetXtra (-fref-1);
- DeleteObject (flaString);
- } else {
- newFormulas << (long)thisFormula;
- references << fref;
- }
- } else {
- isGood = false;
+ if (! process_formula (thisFormula)) {
break;
}
+
} else {
references << -1;
}
@@ -2686,6 +2682,7 @@ void _Matrix::MakeMeSimple (void) {
cmd = new _CompiledMatrixData;
cmd->has_volatile_entries = false;
+ cmd->stackDepth = stackLength;
for (unsigned long k = 0; k < newFormulas.lLength; k++) {
cmd->has_volatile_entries = ((_Formula*)newFormulas.get(k))->ConvertToSimple(varList) || cmd->has_volatile_entries;
@@ -5359,7 +5356,7 @@ void _Matrix::CompressSparseMatrix (bool transpose, hyFloat * stash) {
_Matrix* _Matrix::Exponentiate (hyFloat scale_to, bool check_transition, _Matrix * existing_storage) {
// find the maximal elements of the matrix
-
+
try {
if (!is_square()) {
throw _String ("Exponentiate is not defined for non-square matrices");
@@ -5679,8 +5676,9 @@ _Matrix* _Matrix::Exponentiate (hyFloat scale_to, bool check_transition, _Mat
if (check_transition) {
bool pass = true;
if (result->is_dense()) {
- for (unsigned long r = 0L; r < result->lDim; r += result->vDim) {
+ for (unsigned long r = 0L; r < result->lDim; r += result->vDim + 1) {
if (result->theData[r] > 1.) {
+ //printf ("\nDIAGONAL > 1 (%d, %g)\n", r, result->theData[r]);
pass = false;
break;
}
=====================================
src/core/tree.cpp
=====================================
@@ -2750,14 +2750,7 @@ void _TheTree::ExponentiateMatrices (_List& expNodes, long tc, long catI
}
-
- //if (hasExpForm && (3*(flatTree.countitems() + flatLeaves.countitems()) < (mes_counter-b4)))
- // printf ("%d/%d (%d)\n\n", mes_counter-b4, expNodes.lLength, flatTree.countitems() + flatLeaves.countitems());
- //printf ("%ld %d\n", nodesToDo.lLength, hasExpForm);
- //ObjectToConsole(&isExplicitForm);
-
- unsigned long id;
-
+
_List * computedExponentials = hasExpForm? new _List (matrixQueue.lLength) : nil;
_SimpleList parallel, serial;
@@ -2773,14 +2766,14 @@ void _TheTree::ExponentiateMatrices (_List& expNodes, long tc, long catI
if (parallel.lLength) {
#ifdef _OPENMP
#if _OPENMP>=201511
- #pragma omp parallel for default(shared) private (id) schedule(dynamic, cs) proc_bind(spread) if (nt>1) num_threads (nt)
+ #pragma omp parallel for default(shared) schedule(dynamic, cs) proc_bind(spread) if (nt>1) num_threads (nt)
#else
#if _OPENMP>=200803
- #pragma omp parallel for default(shared) private (id) schedule(guided) proc_bind(spread) if (nt>1) num_threads (nt)
+ #pragma omp parallel for default(shared) schedule(guided) proc_bind(spread) if (nt>1) num_threads (nt)
#endif
#endif
#endif
- for (id = 0; id < parallel.lLength; id++) {
+ for (long id = 0L; id < parallel.lLength; id++) {
long matrixID = parallel.get (id);
if (isExplicitForm.list_data[matrixID] == 0 || !hasExpForm) { // normal matrix to exponentiate
((_CalcNode*) nodesToDo(matrixID))->SetCompExp ((_Matrix*)matrixQueue(matrixID), catID, true);
@@ -2790,27 +2783,12 @@ void _TheTree::ExponentiateMatrices (_List& expNodes, long tc, long catI
}
}
- for ( id = 0; id < serial.lLength; id++) {
+ for (long id = 0L; id < serial.lLength; id++) {
long matrixID = serial.get (id);
_Matrix *already_computed = ((_Matrix*)matrixQueue(matrixID));
(*computedExponentials) [matrixID] = already_computed;
}
- /*for (matrixID = 0; matrixID < matrixQueue.lLength; matrixID++) {
- if (isExplicitForm.list_data[matrixID] == 0 || !hasExpForm) { // normal matrix to exponentiate
- ((_CalcNode*) nodesToDo(matrixID))->SetCompExp ((_Matrix*)matrixQueue(matrixID), catID, true);
- } else {
- if (isExplicitForm.list_data[matrixID] > 0) {
- (*computedExponentials) [matrixID] = ((_Matrix*)matrixQueue(matrixID))->Exponentiate(1., true);
- } else {
- _Matrix *already_computed = ((_Matrix*)matrixQueue(matrixID));
- (*computedExponentials) [matrixID] = already_computed;
- //printf ("\tNO MATRIX UPDATE %d (%d)\n", matrixID, already_computed->GetReferenceCounter());
- //ObjectToConsole( (*computedExponentials) [matrixID] );
- //already_computed ->AddAReference();
- }
- }
- }*/
if (computedExponentials) {
_CalcNode * current_node = nil;
View it on GitLab: https://salsa.debian.org/med-team/hyphy/-/commit/14ff91d5b3828bdefe278130f1e1f97d1472ef2c
--
View it on GitLab: https://salsa.debian.org/med-team/hyphy/-/commit/14ff91d5b3828bdefe278130f1e1f97d1472ef2c
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/20240113/7b665db7/attachment-0001.htm>
More information about the debian-med-commit
mailing list