Test showing error

pounders/py/formquad.py

@@ -9,7 +9,7 @@
 
 def formquad(X, F, delta, xk_in, np_max, Pars, vf):
     """
-    formquad(X, F, delta, xk_in, np_max, Pars, vf) -> [Mdir, np, valid, G, H, Mind]
+    formquad(X, F, delta, xk_in, np_max, Pars, vf) -> [Mdir, mp, valid, G, H, Mind]
     Computes the parameters for m quadratics
         ### FIX COMMENT Line 15 ###
         Q_i(x) = C(i) + G(:,i)'*x + 0.5*x'*H(:,:,i)*x,  i=1:m
@@ -31,7 +31,7 @@ def formquad(X, F, delta, xk_in, np_max, Pars, vf):
     vf      [log] Flag indicating you just want to check model validity
     --OUTPUTS----------------------------------------------------------------
     Mdir    [dbl] [(n-np+1)-by-n]  Unit directions to improve model
-    np      [int] Number of interpolation points (=length(Mind))
+    mp      [int] Number of interpolation points (=length(Mind))
     valid   [log] Flag saying if model is valid within Pars[2]*delta
     G       [dbl] [n-by-m]  Matrix of model gradients centered at X[xkin]
     H       [dbl] [n-by-n-by-m]  Array of model Hessians centered at X[xkin]

pounders/py/pounders.py

@@ -338,13 +338,6 @@ def pounders(Ffun, X_0, n, nf_max, g_tol, delta_0, m, Low, Upp, Prior=None, Opti
             # Need to check because model may be valid after Xsp evaluation
             [Mdir, mp, valid, _, _, _] = formquad(X[: nf + 1, :], F[: nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], 1)
             if not valid:  # ! One strategy for choosing model-improving point:
-                # Update model (exists because delta & xk_in unchanged)
-                D = X[: nf + 1] - X[xk_in]
-                Res[: nf + 1, :] = (F[: nf + 1, :] - Cres) - np.diagonal(0.5 * D @ (np.tensordot(D, Hres, axes=1))).T
-                [_, _, valid, Gres, Hresdel, Mind] = formquad(X[: nf + 1, :], Res[: nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], False)
-                Hres = Hres + Hresdel
-                # Update for modelimp; Cres unchanged b/c xk_in unchanged
-                G, H = combinemodels(Cres, Gres, Hres)
                 # Evaluate model-improving points to pick best one
                 # May eventually want to normalize Mdir first for infty norm
                 # Plus directions
@@ -383,6 +376,13 @@ def pounders(Ffun, X_0, n, nf_max, g_tol, delta_0, m, Low, Upp, Prior=None, Opti
                     # Don't actually use
                     for j in range(m):
                         Gres[:, j] = Gres[:, j] + Hres[:, :, j] @ D.T
+                # Update model (exists because delta & xk_in unchanged)
+                D = X[: nf + 1] - X[xk_in]
+                Res[: nf + 1, :] = (F[: nf + 1, :] - Cres) - np.diagonal(0.5 * D @ (np.tensordot(D, Hres, axes=1))).T
+                [_, _, valid, Gres, Hresdel, Mind] = formquad(X[: nf + 1, :], Res[: nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], False)
+                Hres = Hres + Hresdel
+                # Update for modelimp; Cres unchanged b/c xk_in unchanged
+                G, H = combinemodels(Cres, Gres, Hres)
     if printf:
         print("Number of function evals exceeded")
     flag = ng

pounders/py/tests/test_showing_error_on_row10.py

@@ -0,0 +1,61 @@
+"""
+Unit test of compute function
+"""
+
+import os
+
+import ibcdfo.pounders as pdrs
+import numpy as np
+from calfun import calfun
+from dfoxs import dfoxs
+
+
+if not os.path.exists("benchmark_results"):
+    os.makedirs("benchmark_results")
+
+dfo = np.loadtxt("/home/mmenickelly/BenDFO/data/dfo.dat")
+
+spsolver = 2
+nf_max = 500
+g_tol = 1e-13
+factor = 10
+printf = False
+combinemodels = pdrs.identity_combine
+hfun = lambda F: np.squeeze(F)
+Opts = {"printf": 1, "spsolver": 1, "hfun": hfun, "combinemodels": combinemodels}
+
+for row, (nprob, n, m, factor_power) in enumerate(dfo[10:11]):
+    n = int(n)
+    m = int(m)
+
+    def Ffun(y):
+        out = calfun(y, m, int(nprob), "smooth", 0, num_outs=2)[0]
+        # assert len(out) == m, "Incorrect output dimension"
+        return np.squeeze(out)
+
+    X_0 = dfoxs(n, nprob, int(factor**factor_power))
+    Low = -np.inf * np.ones((1, n))  # 1-by-n Vector of lower bounds [zeros(1, n)]
+    Upp = np.inf * np.ones((1, n))  # 1-by-n Vector of upper bounds [ones(1, n)]
+    nfs = 1
+    F_init = np.zeros((1, 1))
+    F_init[0] = Ffun(X_0)
+    xind = 0
+    delta = 0.1
+
+    Prior = {"X_init": X_0, "F_init": F_init, "nfs": nfs, "xk_in": xind}
+
+    Results = {}
+
+    filename = "./benchmark_results/pounders4py_nf_max=" + str(nf_max) + "_prob=" + str(row) + "_spsolver=" + str(spsolver) + "_hfun=" + combinemodels.__name__ + ".mat"
+    Prior = {"nfs": 1, "F_init": F_init, "X_init": X_0, "xk_in": xind}
+
+    [X, F, hF, flag, xk_best] = pdrs.pounders(Ffun, X_0, n, nf_max, g_tol, delta, 1, Low, Upp, Prior=Prior, Options=Opts, Model={})
+
+    evals = F.shape[0]
+
+    assert flag != 1, "pounders failed"
+    assert hfun(F[0]) > hfun(F[xk_best]), "No improvement found"
+    assert X.shape[0] <= nf_max + nfs, "POUNDERs grew the size of X"
+
+    if flag == 0:
+        assert evals <= nf_max + nfs, "POUNDERs evaluated more than nf_max evaluations"