feat: Improve numerical stability and performance

This commit introduces several improvements to the codebase:

- **Numerical Stability:**
  - Adds a mechanism to enforce constraint violations during likelihood function computations, controlled by the `ENFORCE_CONSTRAINT_VIOLATIONS` environment variable.
  - Improves warnings for numerical issues in `GradientLocateTheBump`.

- **Performance:**
  - Optimizes matrix multiplication functions using SIMD instructions for better performance.

- **Build:**
  - Adds the `-Wno-maybe-uninitialized` compiler flag for GNU compilers.
  - Ensures `HYPHYMPI` is always built and installed if MPI is found.

- **Version:**
  - Bumps the version to 2.5.82.

Author: spond

CMakeLists.txt

@@ -113,6 +113,7 @@ endif()
 
 if(CMAKE_CXX_COMPILER_ID MATCHES "GNU")
     list(APPEND COMMON_COMPILE_FLAGS "-Wno-overloaded-virtual")
+    list(APPEND COMMON_COMPILE_FLAGS "-Wno-maybe-uninitialized")
 endif()
 
 if(CMAKE_CXX_COMPILER_ID MATCHES "GNU" OR CMAKE_CXX_COMPILER_ID MATCHES "Clang")
@@ -255,7 +256,7 @@ install(TARGETS hyphy RUNTIME DESTINATION bin)
 # Other Executables
 #------------------------------------------------------------------------------
 if(MPI_FOUND)
-    add_executable(HYPHYMPI EXCLUDE_FROM_ALL ${SRC_COMMON} ${SRC_UNIXMAIN})
+    add_executable(HYPHYMPI ${SRC_COMMON} ${SRC_UNIXMAIN})
     target_compile_options(HYPHYMPI PRIVATE ${COMMON_COMPILE_FLAGS} ${HYPHY_SIMD_FLAGS})
     target_compile_definitions(HYPHYMPI PRIVATE
         __AFYP_REWRITE_BGM__
@@ -285,7 +286,7 @@ if(MPI_FOUND)
         target_include_directories(HYPHYMPI PRIVATE ${ZLIB_INCLUDE_DIRS})
         target_compile_definitions(HYPHYMPI PRIVATE __ZLIB__)
     endif()
-    install(TARGETS HYPHYMPI RUNTIME DESTINATION bin OPTIONAL)
+    install(TARGETS HYPHYMPI RUNTIME DESTINATION bin)
 endif()
 
 add_executable(HYPHYDEBUG EXCLUDE_FROM_ALL ${SRC_COMMON} ${SRC_UNIXMAIN})

compile_commands.json

@@ -1,4 +1,4 @@
-RequireVersion ("2.5.21");
+RequireVersion ("2.5.82");
 
 LoadFunctionLibrary("libv3/all-terms.bf"); // must be loaded before CF3x4
 
@@ -31,6 +31,7 @@ LoadFunctionLibrary("libv3/models/rate_variation.bf");
 utility.SetEnvVariable ("NORMALIZE_SEQUENCE_NAMES", TRUE);
 utility.SetEnvVariable ("ASSUME_REVERSIBLE_MODELS", TRUE);
 utility.SetEnvVariable ("USE_MEMORY_SAVING_DATA_STRUCTURES", 1e8);
+utility.SetEnvVariable ("ENFORCE_CONSTRAINT_VIOLATIONS", TRUE);
 
 //relax.OPTIMIZATION_LOGS = 1;
 
@@ -482,7 +483,6 @@ if (relax.model_set == "All") { // run all the models
         relax.distribution          = models.codon.BS_REL.ExtractMixtureDistribution(relax.ge.bsrel_model);
         relax.weight_multipliers    = parameters.helper.stick_breaking (utility.SwapKeysAndValues(utility.MatrixToDict(relax.distribution["weights"])),None);
         relax.constrain_parameters   = parameters.ConstrainMeanOfSet(relax.distribution["rates"],relax.weight_multipliers,1,"relax");
-        
 
 
         relax.i = 0;
@@ -1071,6 +1071,9 @@ function relax.FitMainTestPair (prompt) {
         relax.save_fit_path = io.PromptUserForFilePath ("Save RELAX model fit to this file ['/dev/null' to skip]");
     } 
 
+    
+    relax.stashLF = estimators.TakeLFStateSnapshot (relax.alternative_model.fit[terms.likelihood_function]);
+
     io.SpoolLFToPath(relax.alternative_model.fit[terms.likelihood_function], relax.save_fit_path);
 
 	if (relax.numbers_of_tested_groups == 2 && relax.analysis_run_mode != relax.kGroupMode) {
@@ -1148,7 +1151,7 @@ function relax.FitMainTestPair (prompt) {
                     {estimators.GetGlobalMLE (relax.alternative_model.fit.take2,terms.relax.k), relax.alternative_model.fit.take2 [terms.fit.log_likelihood]}
                 });
 
-			if (relax.alternative_model.fit.take2 [terms.fit.log_likelihood] > relax.alternative_model.fit.take2[terms.fit.log_likelihood]) {
+			if (relax.alternative_model.fit.take2 [terms.fit.log_likelihood] > relax.alternative_model.fit [terms.fit.log_likelihood]) {
 			    relax.stash_fitted_K = relax.fitted.K;
 
 				io.ReportProgressMessageMD("RELAX", "alt-2", "\n### Potential for highly unreliable K inference due to multiple local maxima in the likelihood function, treat results with caution ");
@@ -1171,6 +1174,8 @@ function relax.FitMainTestPair (prompt) {
 
 				relax.alternative_model.fit = relax.alternative_model.fit.take2;
                 io.SpoolLFToPath(relax.alternative_model.fit.take2[terms.likelihood_function], relax.save_fit_path);
+                relax.stashLF = estimators.TakeLFStateSnapshot (relax.alternative_model.fit[terms.likelihood_function]);
+
 			} else {
 			    estimators.RestoreLFStateFromSnapshot(relax.alternative_model.fit[terms.likelihood_function], relax.take1_snapshot);
 			}

res/TemplateBatchFiles/SelectionAnalyses/RELAX.bf

@@ -3120,7 +3120,7 @@ bool _ElementaryCommand::HandleSetParameter(_ExecutionList &current_program) {
       _SimpleList _modelCache;
       _AVLListXL modelCache(&_modelCache);
 
-      for (long i = 0; i < node_names.countitems(); i++) {
+      for (unsigned long i = 0; i < node_names.countitems(); i++) {
         const _String object_name = *(_String *)node_names.GetItem(i);
         _CalcNode *tree_node = nil;
 

src/core/batchlanruntime.cpp

@@ -115,7 +115,7 @@ _String const kEmptyString, kPromptForFilePlaceholder("PROMPT_FOR_FILE"),
                     "\"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.81"),
+    kHyPhyVersion = _String("2.5.82"),
 
     kNoneToken = "None", kNullToken = "null",
     kNoKWMatch = "__input_value_not_given__",

src/core/global_things.cpp

@@ -351,9 +351,11 @@ _String const accept_branch_lengths("ACCEPT_BRANCH_LENGTHS"),
     use_traversal_heuristic("USE_TRAVERSAL_HEURISTIC"),
     // TODO (20170413): don't remember what this does; , see @
     // _DataSetFilter::MatchStartNEnd #DEPRECATE
-    verbosity_level_string("VERBOSITY_LEVEL")
+    verbosity_level_string("VERBOSITY_LEVEL"),
     // controls verbosity level during optimization and other long-running
-    // operations
+    // operations,
+    enforce_constraint_violations("ENFORCE_CONSTRAINT_VIOLATIONS")
+
     ;
 
 /** default values get set up in hy_global:: */

src/core/hbl_env.cpp

@@ -49,43 +49,53 @@ namespace hy_env {
  * @brief Check if the value of the environment variable is "true".
  *
  * @param name The name of the environment variable (in HBL).
- * @param default_true The default value to return if the variable is not defined.
- * @return true if the variable is defined and true-like (not zero), false otherwise.
+ * @param default_true The default value to return if the variable is not
+ * defined.
+ * @return true if the variable is defined and true-like (not zero), false
+ * otherwise.
  */
 bool EnvVariableTrue(_String const &name, bool default_true = false);
 
 /**
  * @brief Look up the default numeric value for a given environment variable.
  *
  * @param name The name of the environment variable (in HBL).
- * @return The default value if available and of the correct type, otherwise HY_INVALID_RETURN_VALUE.
+ * @return The default value if available and of the correct type, otherwise
+ * HY_INVALID_RETURN_VALUE.
  */
 hyFloat EnvVariableGetDefaultNumber(_String const &name);
 
 /**
- * @brief Look up the default value for a given environment variable, checking that it is of a particular type.
+ * @brief Look up the default value for a given environment variable, checking
+ * that it is of a particular type.
  *
  * @param name The name of the environment variable (in HBL).
  * @param type The expected type of the variable.
- * @return The default value if available and of the correct type, otherwise nil.
+ * @return The default value if available and of the correct type, otherwise
+ * nil.
  */
 HBLObjectRef EnvVariableGetDefault(_String const &name, unsigned long type);
 
 /**
- * @brief Look up the value for a given environment variable, checking that it is of a particular type.
+ * @brief Look up the value for a given environment variable, checking that it
+ * is of a particular type.
  *
  * @param name The name of the environment variable (in HBL).
  * @param type The expected type of the variable.
- * @return The current value if set and of the correct type, otherwise the default value if available and of the correct type, otherwise nil.
+ * @return The current value if set and of the correct type, otherwise the
+ * default value if available and of the correct type, otherwise nil.
  */
 HBLObjectRef EnvVariableGet(_String const &name, unsigned long type);
 
 /**
  * @brief Look up the numeric value for a given environment variable.
  *
  * @param name The name of the environment variable (in HBL).
- * @param default_value The default value to return if the variable is not defined.
- * @return The current value if set and of the correct type, otherwise the default value if available and of the correct type, otherwise HY_INVALID_RETURN_VALUE.
+ * @param default_value The default value to return if the variable is not
+ * defined.
+ * @return The current value if set and of the correct type, otherwise the
+ * default value if available and of the correct type, otherwise
+ * HY_INVALID_RETURN_VALUE.
  */
 hyFloat EnvVariableGetNumber(_String const &name,
                              hyFloat default_value = HY_INVALID_RETURN_VALUE);
@@ -100,7 +110,8 @@ hyFloat EnvVariableGetNumber(_String const &name,
 void EnvVariableSet(_String const &name, HBLObjectRef value, bool copy);
 
 /**
- * @brief Set the value for an environment variable (in a given namespace if provided).
+ * @brief Set the value for an environment variable (in a given namespace if
+ * provided).
  *
  * @param name The name of the environment variable (in HBL).
  * @param value The value to set the variable to.
@@ -139,8 +150,9 @@ extern const _String accept_rooted_trees, accept_branch_lengths,
     soft_fileio_exceptions, last_fileio_exception, last_raw_file_prompt,
     include_model_spec, lf_convergence_criterion, try_numeric_sequence_match,
     short_mpi_return, kSCFGCorpus, verbosity_level_string,
-    tolerate_numerical_errors, tolerate_constraint_violation, number_threads,
-    gzip_compression_level, tree_parser_namespace, global_kwargs;
+    enforce_constraint_violations, tolerate_numerical_errors,
+    tolerate_constraint_violation, number_threads, gzip_compression_level,
+    tree_parser_namespace, global_kwargs;
 
 ;
 

src/core/include/hbl_env.h

@@ -2030,6 +2030,9 @@ bool _LikelihoodFunction::SendOffToMPI(long index, void *request) {
 //_______________________________________________________________________________________
 
 bool _LikelihoodFunction::PreCompute(void) {
+
+  bool validated_contstraints = true;
+
   if (compiled_constraints) {
     // populate all the independent variables
     for (unsigned long i = 0; i < compiled_constraints->varIndex.lLength; i++) {
@@ -2082,6 +2085,9 @@ bool _LikelihoodFunction::PreCompute(void) {
       if (!cornholio->IsValueInBounds(tp)) {
         ReportWarning(_String("Failing bound checks on ") &
                       *cornholio->GetName() & " = " & _String(tp, "%25.16g"));
+        ReportWarning(_String(
+            cornholio->GetFormulaString(kFormulaStringConversionReportRanges)));
+        validated_contstraints = false;
         // break;
       }
     }
@@ -2099,6 +2105,12 @@ bool _LikelihoodFunction::PreCompute(void) {
     }
     */
 
+    if (!CheckEqual(hy_env::EnvVariableGetNumber(
+                        hy_env::enforce_constraint_violations, 0.),
+                    0.0)) {
+      return validated_contstraints;
+    }
+
     return (i == arrayToCheck->lLength);
   }
 }
@@ -7903,14 +7915,20 @@ hyFloat _LikelihoodFunction::GradientLocateTheBump(hyFloat gPrecision,
       if (maxSoFar < initialValue &&
           !CheckEqual(maxSoFar, initialValue,
                       kMachineEpsilon * errorTolerance)) {
-        _TerminateAndDump(_String(" _LikelihoodFunction::GradientLocateTheBump:"
-                                  " in the Brent loop iteration ") &
-                          long(outcome) & ". " & _String(maxSoFar, "%18.16g") &
-                          " / " & _String(initialValue, "%18.16g") &
-                          ".\n Optimization direction: \n" &
-                          _String((_String *)gradient.toStr()));
-
-        return 0.;
+        _String errorStr =
+            _String(" _LikelihoodFunction::GradientLocateTheBump:"
+                    " in the Brent loop iteration ") &
+            long(outcome) & ". " & _String(maxSoFar, "%18.16g") & " / " &
+            _String(initialValue, "%18.16g") &
+            ".\n Optimization direction: \n" &
+            _String((_String *)gradient.toStr());
+
+        if (hy_env::EnvVariableTrue(hy_env::tolerate_numerical_errors)) {
+          ReportWarningConsole(errorStr);
+        } else {
+          _TerminateAndDump(errorStr);
+          return 0.;
+        }
       }
     }
   }

src/core/likefunc.cpp

@@ -396,7 +396,7 @@ HBLObjectRef _MathObject::ExecuteSingleOp(long opCode, _List *arguments,
     default:
       WarnNotDefined(this, opCode, context);
     }
-  } catch (const _String err) {
+  } catch (const _String &err) {
     if (context) {
       context->ReportError(err);
     } else {

src/core/mathobj.cpp

@@ -5670,7 +5670,7 @@ void _Matrix::CopyMatrixToDenseAndMultiply(_Matrix const &source, hyFloat C,
     } else {
       throw(_String("The target argument must be a dense numeric matrix"));
     }
-  } catch (_String const err) {
+  } catch (_String const &err) {
     HandleApplicationError(_String("Internal error in ") &
                            _String(__PRETTY_FUNCTION__).Enquote() & ": " & err);
   }