Warn about non-idempotent values in MFILE

process/caller.py

@@ -6,6 +6,8 @@
 from process.io.mfile import MFile
 from process.utilities.f2py_string_patch import f2py_compatible_to_string
 from typing import Union, Tuple, TYPE_CHECKING
+import warnings
+from tabulate import tabulate
 
 if TYPE_CHECKING:
     from process.main import Models
@@ -41,7 +43,11 @@ def check_agreement(
         :return: whether values agree or not
         :rtype: bool
         """
-        return np.allclose(previous, current, rtol=1.0e-6)
+        # Check for same shape: mfile length can change between iterations
+        if isinstance(previous, float) or previous.shape == current.shape:
+            return np.allclose(previous, current, rtol=1.0e-6, equal_nan=True)
+        else:
+            return False
 
     def call_models(self, xc: np.ndarray, m: int) -> Tuple[float, np.ndarray]:
         """Evalutate models until results are idempotent.
@@ -113,7 +119,7 @@ def call_models_and_write_output(self, xc: np.ndarray, ifail: int) -> None:
         """
         # TODO The only way to ensure idempotence in all outputs is by comparing
         # mfiles at this stage
-        previous_mfile_arr = None
+        previous_mfile_data = None
 
         try:
             # Evaluate models up to 10 times; any more implies non-converging values
@@ -131,24 +137,36 @@ def call_models_and_write_output(self, xc: np.ndarray, ifail: int) -> None:
                     + "IDEM_MFILE.DAT"
                 )
                 mfile = MFile(mfile_path)
-                mfile_data = {}
-                for var in mfile.data.keys():
-                    mfile_data[var] = mfile.data[var].get_scan(-1)
-
-                # Extract floats from mfile dict into array for straightforward
-                # comparison: only compare floats
-                current_mfile_arr = np.array(
-                    [val for val in mfile_data.values() if isinstance(val, float)]
-                )
-                if previous_mfile_arr is None:
+                # Create mfile dict of float values: only compare floats
+                mfile_data = {
+                    var: val
+                    for var in mfile.data.keys()
+                    if isinstance(val := mfile.data[var].get_scan(-1), float)
+                }
+
+                if previous_mfile_data is None:
                     # First run: need another run to compare with
                     logger.debug(
                         "New mfile created: evaluating models again to check idempotence"
                     )
-                    previous_mfile_arr = np.copy(current_mfile_arr)
+                    previous_mfile_data = mfile_data.copy()
                     continue
 
-                if self.check_agreement(previous_mfile_arr, current_mfile_arr):
+                # Compare previous and current mfiles for agreement
+                nonconverged_vars = {}
+                for var in previous_mfile_data.keys():
+                    previous_value = previous_mfile_data[var]
+                    current_value = mfile_data.get(var, np.nan)
+                    if self.check_agreement(previous_value, current_value):
+                        continue
+                    else:
+                        # Value has changed between previous and current mfiles
+                        nonconverged_vars[var] = [
+                            previous_value,
+                            current_value,
+                        ]
+
+                if len(nonconverged_vars) == 0:
                     # Previous and current mfiles agree (idempotent)
                     logger.debug("Mfiles idempotent, returning")
                     # Divert OUT.DAT and MFILE.DAT output back to original files
@@ -158,17 +176,37 @@ def call_models_and_write_output(self, xc: np.ndarray, ifail: int) -> None:
                     finalise(self.models, ifail)
                     return
 
-                # Mfiles not yet idempotent: re-evaluate models
+                # Mfiles not yet idempotent: need to re-evaluate models
                 logger.debug("Mfiles not idempotent, evaluating models again")
-                previous_mfile_arr = np.copy(current_mfile_arr)
+                previous_mfile_data = mfile_data.copy()
 
-            raise RuntimeError(
+            # Values haven't all stabilised after 10 evaluations
+            # Which variables are still changing?
+            non_idempotent_warning = (
                 "Model evaluations at the current optimisation parameter vector "
                 "don't produce idempotent values in the final output."
             )
+            non_idempotent_table = tabulate(
+                [[k, v[0], v[1]] for k, v in nonconverged_vars.items()],
+                headers=["Variable", "Previous value", "Current value"],
+            )
+
+            warnings.warn(
+                f"\033[93m{non_idempotent_warning}\n{non_idempotent_table}\033[0m"
+            )
+
+            # Close idempotence files, write final output file and mfile
+            ft.init_module.close_idempotence_files()
+            finalise(
+                self.models,
+                ifail,
+                non_idempotent_msg=non_idempotent_warning + "\n" + non_idempotent_table,
+            )
+            return
+
         except Exception:
-            # If exception in model evaluations or idempotence can't be
-            # achieved, delete intermediate idempotence files to clean up
+            # If exception in model evaluations delete intermediate idempotence
+            # files to clean up
             ft.init_module.close_idempotence_files()
             raise
 

process/final.py

@@ -3,9 +3,10 @@
 from process import fortran as ft
 from process.fortran import final_module as fm
 from process import output as op
+from process.fortran import process_output as po
 
 
-def finalise(models, ifail):
+def finalise(models, ifail, non_idempotent_msg: None | str = None):
     """Routine to print out the final point in the scan.
 
     Writes to OUT.DAT and MFILE.DAT.
@@ -14,13 +15,20 @@ def finalise(models, ifail):
     :type models: process.main.Models
     :param ifail: error flag
     :type ifail: int
+    :param non_idempotent_msg: warning about non-idempotent variables, defaults to None
+    :type non_idempotent_msg: None | str, optional
     """
     fm.final_header(ifail)
 
     # Output relevant to no optimisation
     if ft.numerics.ioptimz == -2:
         fm.no_optimisation()
 
+    # Print non-idempotence warning to OUT.DAT only
+    if non_idempotent_msg:
+        po.oheadr(ft.constants.nout, "NON-IDEMPOTENT VARIABLES")
+        po.ocmmnt(ft.constants.nout, non_idempotent_msg)
+
     # Write output to OUT.DAT and MFILE.DAT
     op.write(models, ft.constants.nout)
 

setup.py

@@ -44,6 +44,7 @@
         "CoolProp>=6.4",
         "matplotlib>=2.1.1",
         "seaborn>=0.12.2",
+        "tabulate",
     ],
     "extras_require": {
         "test": ["pytest>=5.4.1", "requests>=2.30", "testbook>=0.4"],