Merge pull request #302 from POptUS/291TestM_Final

Finalize benchmark test script

Author: jared321

pounders/py/tests/TestPoundersExtensive.py

@@ -38,6 +38,11 @@ def test_benchmark_pounders(self):
         factor = 10
 
         for row, (nprob, n, m, factor_power) in enumerate(dfo):
+            # TODO: Set nf_max to match values used in MATLAB to allow for
+            # direct comparison.  I suspect that many optimizations are running
+            # down to delta_min.  Therefore, we don't need a special nf_max,
+            # but can add a check to confirm that at least one optimization did
+            # run down to delta_min.
             if row == 0:
                 nf_max = 500  # Testing delta_min stopping on first problem
             else:
@@ -100,6 +105,7 @@ def Ffun_batch(Y):
                 assert hfun_name.startswith("h_")
                 hfun_name = hfun_name.lstrip("h_")
 
+                # TODO: Need to make problem number 1-based to match filenaming scheme in MATLAB
                 filename = RESULT_PATH.joinpath("pounders_nf_max=" + str(nf_max) + "_prob=" + str(row) + "_spsolver=" + str(spsolver) + "_hfun=" + hfun_name + ".mat")
                 Opts = {"printf": printf, "spsolver": spsolver, "hfun": hfun, "combinemodels": combinemodels}
                 Prior = {"nfs": 1, "F_init": F_init, "X_init": X_0, "xk_in": xind}
@@ -130,6 +136,13 @@ def Ffun_batch(Y):
                 # implementations because, for instance, the MATLAB results
                 # store the best approximation index as 1-based as opposed to
                 # 0-based as this test does.
+                #
+                # TODO: Need to make problem number 1-based to match MATLAB's
+                # value.  Normally there should be no need to adjust this since
+                # it's an internal value and we could adjust as needed when
+                # loading the data when alg == POUNDERS_Py.  However, we are
+                # forced to use a 1-based problem number in the filename, so we
+                # should make the value here match the value in the filename.
                 Results = {"alg": "POUNDERS_Py", "problem": "problem " + str(row) + " from More/Wild", "Fvec": F, "H": hF, "X": X, "flag": flag, "xk_best": xk_best}
                 # oct2py.kill_octave() # This is necessary to restart the octave instance,
                 #                      # and thereby remove some caching of inside of oct2py,

pounders/py/tests/compare_results.py

@@ -1,83 +1,29 @@
-import scipy.io
-
 import numpy as np
 
 from .load_results import load_results
 
+_FLAG_DELTA_MIN = -6
 
-def _load_results_m_v2(filename):
-    """
-    POUNDERS/MATLAB v2 format established at commit 4336b866.
-    """
-    EXPECTED_KEYS = {"alg", "problem", "H", "Fvec", "X", "flag", "xk_best"}
-
-    contents = scipy.io.loadmat(filename)
-    keys = [k for k in contents.keys() if not k.startswith("__")]
-    assert len(keys) == 1
-    assert keys[0] == "Results"
-
-    # Only one valid set of results across all three known hfun cases.
-    data = None
-    tmp = contents[keys[0]]
-    assert len(tmp) == 3
-    for hfun in range(len(tmp)):
-        # See if we can find that one valid result for this hfun case.
-        tmp_i = [e for e in tmp[hfun] if np.squeeze(e).ndim == 0]
-        if tmp_i:
-            assert len(tmp_i) == 1
-            assert data is None
-            data = tmp_i[0][0]
-            assert data is not None
-
-    assert set(data.dtype.names) == EXPECTED_KEYS
-
-    algorithm = data["alg"][0][0]
-    problem = data["problem"][0][0]
-
-    H = np.squeeze(data["H"][0])
-    assert H.ndim == 1
-    n_evaluations = len(H)
-    assert all(np.isreal(H))
-
-    Fvec = np.squeeze(data["Fvec"][0])
-    assert Fvec.ndim == 2
-    tmp, _ = Fvec.shape
-    assert tmp == n_evaluations
-    assert all(np.isreal(Fvec.flatten()))
-
-    X = np.squeeze(data["X"][0])
-    assert X.ndim == 2
-    tmp, _ = X.shape
-    assert tmp == n_evaluations
-    assert all(np.isreal(X.flatten()))
-    assert all(np.isfinite(X.flatten()))
-
-    flag = np.squeeze(data["flag"][0])
-    assert np.isreal(flag)
-    assert np.isfinite(flag)
-    # Stored as 1-based index, but needs to be 0-based index for working with
-    # Python arrays in this code.
-    xk_best = np.squeeze(data["xk_best"][0]) - 1
-    assert xk_best in range(0, len(H))
-
-    return algorithm, problem, X, Fvec, H, xk_best, flag
+
+def _failed(flag):
+    # Having the optimization terminate due to the trust region radius
+    # shrinking down to delta_min does not necessarily indicate a failure.  If
+    # delta_min is well-specified, it could be treated as a success.
+    return (flag < 0) and (flag != _FLAG_DELTA_MIN)
 
 
 def compare_results(filename_benchmark, filename_result):
     """
     .. todo::
         * Allow for users to specify nonzero tolerances if the use case arises.
-        * Allow for checking Python and MATLAB results on a set of problems on
-          which we expect all optimizations to find the same local minimizer.
-          This would require nonzero tolerances.
 
     :param filename_benchmark: Filename of |pounders| ``.mat``-format
         benchmarking result that calling code considers to be the accepted
         reference.
     :param filename_result: Filename of |pounders| ``.mat``-format benchmarking
         result that calling code wishes to check against the reference.
     :return: True if the files correspond to identical test setups and contain
-        bitwise-identical results.
+        valid, bitwise-identical results.
     """
     # ----- HARDCODED VALUES
     RED = "\033[0;91;1m"  # Bright Red/bold
@@ -95,7 +41,7 @@ def error(msg):
         error(f"New result has different filename ({filename_result.stem})")
         return False
 
-    ref_alg, ref_problem, X_ref, F_ref, H_ref, x_best_ref, flag_ref = _load_results_m_v2(filename_benchmark)
+    ref_alg, ref_problem, X_ref, F_ref, H_ref, x_best_ref, flag_ref = load_results(filename_benchmark)
     if not all(np.isfinite(H_ref)):
         error("Non-finite h values in benchmark")
         return False
@@ -111,37 +57,27 @@ def error(msg):
         error("Non-finite Fvec values in new results")
         return False
 
-    # TODO: Once we are testing v3 against v3, remove this check since
-    # load_result() error checks the algorithm name and we would like to check
-    # MATLAB results against Python results.
-    if ref_alg not in ["POUNDERs"]:
-        error(f"Invalid algorithm name ({ref_alg}) for benchmark")
-        return False
-    elif new_alg != "POUNDERS_M":
-        msg = "Benchmark and new result used different algorithms ({} != {})"
-        error(msg.format(ref_alg, new_alg))
-        return False
-
     if new_problem != ref_problem:
         msg = "Benchmark and new result solve different problems ({} != {})"
         error(msg.format(ref_problem, new_problem))
         return False
 
     # ----- COMPARE NEW RESULTS AGAINST BENCHMARK
-    if len(H_new) != len(H_ref):
-        error(f"H arrays have different lengths ({len(H_ref)} != {len(H_new)})")
-        return False
-    assert F_new.shape == F_ref.shape
-    assert X_new.shape == X_ref.shape
-
-    # Don't fail immediately if values are different so that we can provide
-    # users with all such differences in one go.
-    msgs = []
+    # These checks are designed under the assumption that the prime use of this
+    # function is to detect if two results are not *identical*.
+    #
+    # Even so, we don't fail immediately if values are different so that we can
+    # provide users with all such differences in one go.
+    assert F_new.shape[1] == F_ref.shape[1]
+    assert X_new.shape[1] == X_ref.shape[1]
+
+    errors = []
+    warnings = []
     if x_best_new != x_best_ref:
-        msgs += [f"Best approximation indices differ ({x_best_new} != {x_best_ref})"]
+        errors += [f"Best approximation indices differ ({x_best_new} != {x_best_ref})"]
     if flag_new != flag_ref:
-        msgs += [f"Flags differ ({flag_new} != {flag_ref})"]
-    if (flag_new >= 0) and (flag_ref >= 0):
+        errors += [f"Flags differ ({flag_new} != {flag_ref})"]
+    if (not _failed(flag_new)) and (not _failed(flag_ref)):
         # Only show comparison if both ran without a hard failure.  For
         # instance, I would like to see the these comparisons if one or both
         # were simply nonconvergent.
@@ -153,19 +89,34 @@ def error(msg):
         F_best_new = F_new[x_best_new]
         H_best_new = H_new[x_best_new]
 
+        if flag_ref == _FLAG_DELTA_MIN:
+            warnings += ["Benchmark reached delta_min"]
+        if flag_new == _FLAG_DELTA_MIN:
+            warnings += ["New result reached delta_min"]
+
         if H_best_new != H_best_ref:
             abs_diff = np.fabs(H_best_new - H_best_ref)
-            msgs += [f"H absolute difference = {abs_diff}"]
+            errors += [f"H absolute difference = {abs_diff}"]
         if any(F_best_new != F_best_ref):
             max_abs_diff = np.max(np.fabs(F_best_new - F_best_ref))
-            msgs += [f"Fvec max absolute difference = {max_abs_diff}"]
+            errors += [f"Fvec max absolute difference = {max_abs_diff}"]
         if any(X_best_new != X_best_ref):
             max_abs_diff = np.max(np.fabs(X_best_new - X_best_ref))
-            msgs += [f"X max absolute difference = {max_abs_diff}"]
-
-    if msgs:
-        error("\n\t".join(msgs))
+            errors += [f"X max absolute difference = {max_abs_diff}"]
+    else:
+        # We've already reported an error if the flags differ and identical
+        # "bad" flags is not necessarily a failure.
+        if _failed(flag_ref):
+            warnings += [f"Benchmark failed with flag={flag_ref}"]
+        if _failed(flag_new):
+            warnings += [f"New result failed with flag={flag_new}"]
+
+    if errors:
+        error("\n\t".join(errors + warnings))
         return False
+    elif warnings:
+        print(f"{BLUE}PASS{NC}\n\t" + "\n\t".join(warnings))
+        return True
 
     print(f"{BLUE}PASS{NC}")
     return True

pounders/py/tests/load_results.py

@@ -1,3 +1,5 @@
+import copy
+
 import scipy.io
 
 import numpy as np
@@ -11,18 +13,27 @@ def load_results(filename):
     """
     EXPECTED_KEYS = {"alg", "problem", "H", "Fvec", "X", "flag", "xk_best"}
 
-    # When adapting this to loading from MATLAB-generated files on GCE, I was
-    # forced to add the simplify_cells argument.  Without this, some loaded
-    # results had a "W" or a "V" field instead of the "X" field.  When I loaded
-    # those same files manually in Python using the exact same venv and
-    # scipy.io.loadmat(filename), I always saw the "X" field and didn't see the
-    # others.  I don't understand why this was happening nor why the addition
-    # of the flag solves it.
-    data = scipy.io.loadmat(filename, simplify_cells=True)
+    with open(filename, "rb") as fptr:
+        # If I call this function twice in a row on the same MATLAB-generated
+        # file (using both scipy v1.16.2 and v1.17.1), on some occasions the "X"
+        # field will appear with a different name (e.g., "W") on the second
+        # call.  This is apparently related to the change to xk_best below.  If
+        # I rerun that same test after exiting Python, I get the same result,
+        # which implies that the file itself is not being altered.  I hope that
+        # this is the case since we load in readonly mode.  This even happens if
+        # the two files loaded back-to-back are different.
+        #
+        # Therefore, I suspect that there is some sort of caching of or
+        # accidental referencing to loaded data.  Caching does not make sense
+        # when loading two different files.  However, using a deepcopy resolves
+        # the issue.  A simple copy is not sufficient, which makes sense since
+        # data is a nested dictionary.
+        data = copy.deepcopy(scipy.io.loadmat(fptr))
     keys = [k for k in data.keys() if not k.startswith("__")]
     assert set(keys) == EXPECTED_KEYS
 
     algorithm = str(np.squeeze(data["alg"]))
+    assert algorithm in ["POUNDERS_M", "POUNDERS_Py"]
 
     problem = str(np.squeeze(data["problem"]))
     if (not problem.startswith("problem")) or (not problem.endswith("from More/Wild")):
@@ -32,37 +43,39 @@ def load_results(filename):
     except Exception:
         raise ValueError(f"Invalid problem spec ({problem})")
 
-    H = np.squeeze(data["H"])
+    H = np.atleast_1d(np.squeeze(data["H"]))
     assert H.ndim == 1
     n_evaluations = len(H)
     assert all(np.isreal(H))
+    # Non-finite can happen in our tests, so checking finiteness of H must be
+    # handled by calling code.
 
     # Fvec could be a scalar at each evaluation
-    Fvec = np.atleast_2d(np.squeeze(data["Fvec"]))
+    Fvec = np.atleast_2d(data["Fvec"])
     assert Fvec.ndim == 2
     tmp, _ = Fvec.shape
     assert tmp == n_evaluations
     assert all(np.isreal(Fvec.flatten()))
+    # Non-finite can happen in our tests, so checking finiteness of Fvec must
+    # be handled by calling code.
 
     # X could be a scalar at each evaluation
-    X = np.atleast_2d(np.squeeze(data["X"]))
+    X = np.atleast_2d(data["X"])
     assert X.ndim == 2
     tmp, _ = X.shape
     assert tmp == n_evaluations
     assert all(np.isreal(X.flatten()))
     assert all(np.isfinite(X.flatten()))
 
     flag = np.squeeze(data["flag"])
-    assert np.isreal(flag)
-    assert np.isfinite(flag)
+    assert (flag >= 0.0) or (flag in [-6, -5, -4, -3, -2, -1])
+
     xk_best = np.squeeze(data["xk_best"])
     if algorithm == "POUNDERS_M":
         # The MATLAB implementation's test suite saves the index of the best
         # approximation as a 1-based integer.  However, we need to adjust it to
         # 0-based since we are returning Python arrays.
         xk_best -= 1
-    elif algorithm != "POUNDERS_Py":
-        raise ValueError(f"Unknown POUNDERS test algorithm {algorithm}")
-    assert xk_best in range(0, len(H))
+    assert xk_best in range(0, n_evaluations)
 
     return algorithm, problem, X, Fvec, H, xk_best, flag