Repository: POptUS/IBCDFO (language: python)

Issue:
removing duplicate debugging calculation
These were accidentally pulled into `main` in #29.


Discussion / hints:

Fix patch:
diff --git a/pounders/m/checkinputss.m b/pounders/m/checkinputss.m
index 7f3c8f3a..5043af29 100644
--- a/pounders/m/checkinputss.m
+++ b/pounders/m/checkinputss.m
@@ -69,7 +69,7 @@
     flag = 0;
 end
 
-% Check vector of initial function values
+% Check matrix of initial function values
 [nfs2, m2] = size(F0);
 if nfs2 < nfs
     disp('  Error: fewer than nfs function values in F0');
@@ -83,6 +83,11 @@
     disp('  Warning: number of starting f values nfs does not match input F0');
     flag = 0;
 end
+if any(any(isnan(F0)))
+    disp("  Error: F0 contains a NaN.");
+    flag = -1;
+    return
+end
 
 % Check starting point
 if xkin > max(nfs, 1) || xkin < 1 || mod(xkin, 1) ~= 0
diff --git a/pounders/m/pounders.m b/pounders/m/pounders.m
index 333e3bc0..7cb306a7 100644
--- a/pounders/m/pounders.m
+++ b/pounders/m/pounders.m
@@ -57,6 +57,7 @@
 %               > 0 exceeded nfmax evals,   flag = norm of grad at final X
 %               = -1 if input was fatally incorrect (error message shown)
 %               = -2 model failure
+%               = -3 if NaN was encountered
 % xkin    [int] Index of point in X representing approximate minimizer
 %
 % --DEPENDS ON-------------------------------------------------------------
@@ -129,6 +130,13 @@
     F = zeros(nfmax, m); % Stores the function values
     nf = 1;
     F(nf, :) = fun(X(nf, :));
+    if any(isnan(F(nf, :)))
+        disp("A NaN was encountered in an objective evaluation. Exiting.");
+        X = X(1:nf, :);
+        F = F(1:nf, :);
+        flag = -3;
+        return
+    end
     if printf
         fprintf('%4i    Initial point  %11.5e\n', nf, sum(F(nf, :).^2));
     end
@@ -144,7 +152,6 @@
 end
 
 Res = zeros(size(F)); % Stores the residuals for model updates
-ResJ = zeros(size(F)); % Stores the residuals for model updates
 Cres = F(xkin, :);
 Hres = zeros(n, n, m);
 % H = zeros(n); G = zeros(n,1); c = Fs(xkin);
@@ -156,15 +163,7 @@
     % 1a. Compute the interpolation set.
     for i = 1:nf
         D = X(i, :) - X(xkin, :);
-        for j = 1:m
-            Res(i, j) = (F(i, j) - Cres(j)) - .5 * D * Hres(:, :, j) * D';
-        end
-        ResJ(i, :) = F(i, :) - Cres - .5 * D * reshape(D * reshape(Hres, n, m * n), n, m);
-        if norm(Res(i, :)) > 1e-16
-            assert(norm(ResJ(i, :) - Res(i, :)) / norm(Res(i, :)) < 1e-12, "A");
-        else
-            assert(norm(ResJ(i, :) - Res(i, :)) < 1e-12, "A");
-        end
+        Res(i, :) = F(i, :) - Cres - .5 * D * reshape(D * reshape(Hres, n, m * n), n, m);
     end
     [Mdir, np, valid, Gres, Hresdel, Mind] = ...
         formquad(X(1:nf, :), Res(1:nf, :), delta, xkin, npmax, Par, 0);
@@ -174,20 +173,19 @@
             nf = nf + 1;
             X(nf, :) = min(U, max(L, X(xkin, :) + Mdir(i, :))); % Temp safeguard
             F(nf, :) = fun(X(nf, :));
-            Fs(nf) = hfun(F(nf, :));
+            if any(isnan(F(nf, :)))
+                disp("A NaN was encountered in an objective evaluation. Exiting.");
+                X = X(1:nf, :);
+                F = F(1:nf, :);
+                flag = -3;
+                return
+            end
+                Fs(nf) = hfun(F(nf, :));
             if printf
                 fprintf('%4i   Geometry point  %11.5e\n', nf, Fs(nf));
             end
             D = Mdir(i, :);
-            for j = 1:m
-                Res(nf, j) = (F(nf, j) - Cres(j)) - .5 * D * Hres(:, :, j) * D';
-            end
-            ResJ(nf, :) = F(nf, :) - Cres - .5 * D * reshape(D * reshape(Hres, n, m * n), n, m);
-            if norm(Res(i, :)) > 1e-16
-                assert(norm(ResJ(nf, :) - Res(nf, :)) / norm(Res(i, :)) < 1e-12, "A");
-            else
-                assert(norm(ResJ(nf, :) - Res(nf, :)) < 1e-12, "A");
-            end
+            Res(nf, :) = F(nf, :) - Cres - .5 * D * reshape(D * reshape(Hres, n, m * n), n, m);
         end
         if nf >= nfmax
             break
@@ -231,6 +229,13 @@
                 nf = nf + 1;
                 X(nf, :) = min(U, max(L, X(xkin, :) + Mdir(i, :))); % Temp safeg.
                 F(nf, :) = fun(X(nf, :));
+                if any(isnan(F(nf, :)))
+                    disp("A NaN was encountered in an objective evaluation. Exiting.");
+                    X = X(1:nf, :);
+                    F = F(1:nf, :);
+                    flag = -3;
+                    return
+                end
                 Fs(nf) = hfun(F(nf, :));
                 if printf
                     fprintf('%4i   Critical point  %11.5e\n', nf, Fs(nf));
@@ -305,6 +310,13 @@
         nf = nf + 1;
         X(nf, :) = Xsp;
         F(nf, :) = fun(X(nf, :));
+        if any(isnan(F(nf, :)))
+            disp("A NaN was encountered in an objective evaluation. Exiting.");
+            X = X(1:nf, :);
+            F = F(1:nf, :);
+            flag = -3;
+            return
+        end
         Fs(nf) = hfun(F(nf, :));
 
         if mdec ~= 0
@@ -350,15 +362,7 @@
             % Update model (exists because delta & xkin unchanged)
             for i = 1:nf
                 D = (X(i, :) - X(xkin, :));
-                for j = 1:m
-                    Res(i, j) = (F(i, j) - Cres(j)) - .5 * D * Hres(:, :, j) * D';
-                end
-                ResJ(i, :) = F(i, :) - Cres - .5 * D * reshape(D * reshape(Hres, n, m * n), n, m);
-                if norm(Res(i, :)) > 1e-16
-                    assert(norm(ResJ(i, :) - Res(i, :)) / norm(Res(i, :)) < 1e-12, "A");
-                else
-                    assert(norm(ResJ(i, :) - Res(i, :)) < 1e-12, "A");
-                end
+                Res(i, :) = F(i, :) - Cres - .5 * D * reshape(D * reshape(Hres, n, m * n), n, m);
             end
             [~, ~, valid, Gres, Hresdel, Mind] = ...
                 formquad(X(1:nf, :), Res(1:nf, :), delta, xkin, npmax, Par, 0);
@@ -390,6 +394,13 @@
             nf = nf + 1;
             X(nf, :) = min(U, max(L, X(xkin, :) + Xsp)); % Temp safeguard
             F(nf, :) = fun(X(nf, :));
+            if any(isnan(F(nf, :)))
+                disp("A NaN was encountered in an objective evaluation. Exiting.");
+                X = X(1:nf, :);
+                F = F(1:nf, :);
+                flag = -3;
+                return
+            end
             Fs(nf) = hfun(F(nf, :));
             if printf
                 fprintf('%4i   Model point     %11.5e\n', nf, Fs(nf));


Test patch:
diff --git a/pounders/m/tests/failing_objective.m b/pounders/m/tests/failing_objective.m
new file mode 100644
index 00000000..96a44ae5
--- /dev/null
+++ b/pounders/m/tests/failing_objective.m
@@ -0,0 +1,9 @@
+function [fvec] = failing_objective(x)
+% This function will produce a NaN in one component on average in 1 out of 10
+% queries, otherwise, it produces returns x.
+
+fvec = x;
+
+if rand < 0.1
+    fvec(1) = NaN;
+end
diff --git a/pounders/m/tests/test_failing_objective.m b/pounders/m/tests/test_failing_objective.m
new file mode 100644
index 00000000..6a257010
--- /dev/null
+++ b/pounders/m/tests/test_failing_objective.m
@@ -0,0 +1,31 @@
+% This tests the ability of the MATLAB pounders implementation to
+% gracefully terminate when receiving a NaN as a component in the vector of
+% objective output.
+
+function [] = test_failing_objective()
+
+spsolver = 1;
+
+addpath('../');
+
+nfmax = 1000;
+gtol = 1e-13;
+n = 3;
+m = 3;
+
+xs = [10; 20; 30];
+npmax = 2 * n + 1;
+L = -Inf(1, n);
+U = Inf(1, n);
+nfs = 0;
+X0 = xs';
+F0 = [];
+xkin = 1;
+delta = 0.1;
+printf = 0;
+
+rand('seed', 1);
+objective = @(x)failing_objective(x);
+
+[X, F, flag, xk_best] = pounders(objective, X0, n, npmax, nfmax, gtol, delta, nfs, m, F0, xkin, L, U, printf, spsolver);
+assert(flag == -3, "No NaN was encountered in this test, but (with high probability) should have been.");