diff --git a/CMakeLists.txt b/CMakeLists.txt index f6ef69c62..53b09a36c 100644 --- a/CMakeLists.txt +++ b/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 diff --git a/res/TemplateBatchFiles/SelectionAnalyses/FADE.bf b/res/TemplateBatchFiles/SelectionAnalyses/FADE.bf index d99b9939a..a7822064e 100644 --- a/res/TemplateBatchFiles/SelectionAnalyses/FADE.bf +++ b/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", diff --git a/res/TemplateBatchFiles/SelectionAnalyses/MEME.bf b/res/TemplateBatchFiles/SelectionAnalyses/MEME.bf index 556311c7c..c73a0680a 100644 --- a/res/TemplateBatchFiles/SelectionAnalyses/MEME.bf +++ b/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); diff --git a/res/TemplateBatchFiles/SelectionAnalyses/error-filter.bf b/res/TemplateBatchFiles/SelectionAnalyses/error-filter.bf index b72e35469..6225f6323 100644 --- a/res/TemplateBatchFiles/SelectionAnalyses/error-filter.bf +++ b/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); } diff --git a/src/core/formula.cpp b/src/core/formula.cpp index 489f8db52..721eb74a0 100644 --- a/src/core/formula.cpp +++ b/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; itheNumber) { - 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; itheNumber->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; + } + } + } - } - } } } diff --git a/src/core/global_things.cpp b/src/core/global_things.cpp index 22a9d3046..4f2365158 100644 --- a/src/core/global_things.cpp +++ b/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", diff --git a/src/core/include/formula.h b/src/core/include/formula.h index b5a994ea6..56c7821d8 100644 --- a/src/core/include/formula.h +++ b/src/core/include/formula.h @@ -80,7 +80,7 @@ class _Formula { _List* resultCache; _Stack theStack; _List theFormula; - long* simpleExpressionStatus; + long* simpleExpressionStatus; /** SLKP: 20200924 Added this shorthand to improve memory locality and speed-up SimpleCompute performance diff --git a/src/core/include/matrix.h b/src/core/include/matrix.h index c4a6e9b4e..c2c677bbe 100644 --- a/src/core/include/matrix.h +++ b/src/core/include/matrix.h @@ -75,6 +75,7 @@ struct _CompiledMatrixData { hyFloat * formulaValues; long * formulaRefs; + long stackDepth; bool has_volatile_entries; _SimpleList varIndex, diff --git a/src/core/include/simplelist.h b/src/core/include/simplelist.h index f1b492784..7c174a0cf 100644 --- a/src/core/include/simplelist.h +++ b/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 diff --git a/src/core/likefunc.cpp b/src/core/likefunc.cpp index 6fa1571a3..19f674a99 100644 --- a/src/core/likefunc.cpp +++ b/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; diff --git a/src/core/matrix.cpp b/src/core/matrix.cpp index f78744818..64f253b7d 100644 --- a/src/core/matrix.cpp +++ b/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; iAmISimple(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; iIsEmpty())) { - 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; } diff --git a/src/core/tree.cpp b/src/core/tree.cpp index c81fd4e69..56504acaa 100644 --- a/src/core/tree.cpp +++ b/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;