Merge pull request #79 from POptUS/bugfix/python_mpmax

Bugfix/python mpmax

pounders/py/checkinputss.py

@@ -1,9 +1,9 @@
 import numpy as np
 
 
-def checkinputss(fun, X0, n, mpmax, nfmax, gtol, delta, nfs, m, F0, xkin, L, U):
+def checkinputss(fun, X0, n, npmax, nfmax, gtol, delta, nfs, m, F0, xkin, L, U):
     """
-    checkinputss(fun,X0,n,mpmax,nfmax,gtol,delta,nfs,m,F0,xkin,L,U) -> [flag,X0,mpmax,F0,L,U]
+    checkinputss(fun,X0,n,npmax,nfmax,gtol,delta,nfs,m,F0,xkin,L,U) -> [flag,X0,npmax,F0,L,U]
     Checks the inputs provided to pounders.
     A warning message is produced if a nonfatal input is given (and the input is changed accordingly).
     An error message (flag=-1) is produced if the pounders cannot continue.
@@ -18,7 +18,7 @@ def checkinputss(fun, X0, n, mpmax, nfmax, gtol, delta, nfs, m, F0, xkin, L, U):
     if not callable(fun):
         print("Error: fun is not a function handle")
         flag = -1
-        return [flag, X0, mpmax, F0, L, U]
+        return [flag, X0, npmax, F0, L, U]
     # Verify X0 is the appropriate size
     X0 = np.atleast_2d(X0)
     F0 = np.atleast_2d(F0)
@@ -35,21 +35,21 @@ def checkinputss(fun, X0, n, mpmax, nfmax, gtol, delta, nfs, m, F0, xkin, L, U):
         else:
             print("Error: np.shape(X0)[1] != n")
             flag = -1
-            return [flag, X0, mpmax, F0, L, U, xkin]
+            return [flag, X0, npmax, F0, L, U, xkin]
     # Check max number of interpolation points
-    if mpmax < n + 1 or mpmax > int(0.5 * (n + 1) * (n + 2)):
-        mpmax = max(n + 1, min(mpmax, int(0.5 * (n + 1) * (n + 2))))
-        print(f"Warning: mpmax not in [n+1, 0.5 * (n+1) * (n+2) using {mpmax}")
+    if npmax < n + 1 or npmax > int(0.5 * (n + 1) * (n + 2)):
+        npmax = max(n + 1, min(npmax, int(0.5 * (n + 1) * (n + 2))))
+        print(f"Warning: npmax not in [n+1, 0.5 * (n+1) * (n+2) using {npmax}")
         flag = 0
     # Check standard positive quantities
     if nfmax < 1:
         print("Error: max number of evaluations is less than 1")
         flag = -1
-        return [flag, X0, mpmax, F0, L, U, xkin]
+        return [flag, X0, npmax, F0, L, U, xkin]
     elif gtol <= 0:
         print(" Error: gtol must be positive")
         flag = -1
-        return [flag, X0, mpmax, F0, L, U, xkin]
+        return [flag, X0, npmax, F0, L, U, xkin]
     elif delta <= 0:
         print("Error: delta must be positive")
     # Check number of starting points
@@ -63,31 +63,31 @@ def checkinputss(fun, X0, n, mpmax, nfmax, gtol, delta, nfs, m, F0, xkin, L, U):
         if nfs2 < nfs:
             print("Error: fewer than nfs function values in F0")
             flag = -1
-            return [flag, X0, mpmax, F0, L, U, xkin]
+            return [flag, X0, npmax, F0, L, U, xkin]
         elif nfs > 1 and m != m2:
             print("Error: F0 does not contain the right number of residuals")
             flag = -1
-            return [flag, X0, mpmax, F0, L, U, xkin]
+            return [flag, X0, npmax, F0, L, U, xkin]
         elif nfs2 > nfs:
             print("Warning: number of starting f values nfs does not match input F0")
             flag = 0
         if np.any(np.isnan(F0)):
             print("Error: F0 contains a NaN.")
             flag = -1
-            return [flag, X0, mpmax, F0, L, U, xkin]
+            return [flag, X0, npmax, F0, L, U, xkin]
 
     # Check starting point
     if (xkin > max(nfs - 1, 0)) or (xkin < 0) or (xkin % 1 != 0):  # FixMe: Check what xkin needs to be...
         print("Error: starting point index not an integer between 0 and nfs-1")
         flag = -1
-        return [flag, X0, mpmax, F0, L, U, xkin]
+        return [flag, X0, npmax, F0, L, U, xkin]
     # Check the bounds
     [nfs2, n2] = np.shape(L)
     [nfs3, n3] = np.shape(U)
     if (n3 != n2) or (nfs2 != nfs3):
         print("Error: bound dimensions inconsistent")
         flag = -1
-        return [flag, X0, mpmax, F0, L, U, xkin]
+        return [flag, X0, npmax, F0, L, U, xkin]
     elif n2 != n and (n2 == 1 and nfs2 == n):
         L = L.T
         U = U.T
@@ -96,15 +96,15 @@ def checkinputss(fun, X0, n, mpmax, nfmax, gtol, delta, nfs, m, F0, xkin, L, U):
     elif n2 != n or nfs2 != 1:
         print("Error: bounds are not 1-by-n vectors")
         flag = -1
-        return [flag, X0, mpmax, F0, L, U, xkin]
+        return [flag, X0, npmax, F0, L, U, xkin]
     if np.min(U - L) <= 0:
         print("Error: must have U > L")
         flag = -1
-        return [flag, X0, mpmax, F0, L, U, xkin]
+        return [flag, X0, npmax, F0, L, U, xkin]
     if np.min([np.min(X0[xkin, :] - L), np.min(U - X0[xkin, :])]) < 0:
         print("Error: starting point outside of bounds (L,U)")
         flag = -1
-        return [flag, X0, mpmax, F0, L, U, xkin]
+        return [flag, X0, npmax, F0, L, U, xkin]
     U = U.squeeze()
     L = L.squeeze()
-    return [flag, X0, mpmax, F0, L, U, xkin]
+    return [flag, X0, npmax, F0, L, U, xkin]

pounders/py/formquad.py

@@ -7,7 +7,7 @@
 # from .flipSignQ import flipSignQ
 
 
-def formquad(X, F, delta, xkin, mpmax, Pars, vf):
+def formquad(X, F, delta, xkin, npmax, Pars, vf):
     """
     formquad(X,F,delta,xkin,npmax,Pars,vf) -> [Mdir,np,valid,G,H,Mind]
     Computes the parameters for m quadratics
@@ -48,7 +48,7 @@ def formquad(X, F, delta, xkin, mpmax, Pars, vf):
     # Precompute the scaled displacements (could be expensive for larger nfmax)
     D = np.zeros((nf, n))  # Scaled displacements
 
-    assert isinstance(mpmax, int), "Must be an integer"
+    assert isinstance(npmax, int), "Must be an integer"
     assert isinstance(xkin, int), "Must be an integer"
 
     D = (X[:nf] - X[xkin]) / delta
@@ -103,9 +103,9 @@ def formquad(X, F, delta, xkin, mpmax, Pars, vf):
     M = np.vstack((np.ones(n + 1), D[Mind].T))
     [Q, R] = np.linalg.qr(M.T, mode="complete")
     # [Q, R] = flipSignQ(Q, R, 0, np.shape(Q)[1]-1)
-    # Now we add points until we have mpmax starting with the most recent ones
+    # Now we add points until we have npmax starting with the most recent ones
     i = nf - 1
-    while mp < mpmax or mpmax == n + 1:
+    while mp < npmax or npmax == n + 1:
         if Nd[i] <= Pars[1] and i not in Mind:
             Ny = np.hstack((N, phi2eval(D[[i], :]).T))
             # Update QR

pounders/py/pounders.py

@@ -9,11 +9,11 @@
 from .prepare_outputs_before_return import prepare_outputs_before_return
 
 
-def pounders(fun, X0, n, mpmax, nfmax, gtol, delta, nfs, m, F0, xkin, L, U, printf=0, spsolver=2, hfun=None, combinemodels=None):
+def pounders(fun, X0, n, npmax, nfmax, gtol, delta, nfs, m, F0, xkin, L, U, printf=0, spsolver=2, hfun=None, combinemodels=None):
     """
     POUNDERS: Practical Optimization Using No Derivatives for sums of Squares
       [X,F,flag,xkin] = ...
-           pounders(fun,X0,n,mpmax,nfmax,gtol,delta,nfs,m,F0,xkin,L,U,printf)
+           pounders(fun,X0,n,npmax,nfmax,gtol,delta,nfs,m,F0,xkin,L,U,printf)
 
     This code minimizes output from a structured blackbox function, solving
     min { f(X)=sum_(i=1:m) F_i(x)^2, such that L_j <= X_j <= U_j, j=1,...,n }
@@ -36,7 +36,7 @@ def pounders(fun, X0, n, mpmax, nfmax, gtol, delta, nfs, m, F0, xkin, L, U, prin
     fun     [f h] Function handle so that fun(x) evaluates F (@calfun)
     X0      [dbl] [max(nfs,1)-by-n] Set of initial points  (zeros(1,n))
     n       [int] Dimension (number of continuous variables)
-    mpmax   [int] Maximum number of interpolation points (>n+1) (2*n+1)
+    npmax   [int] Maximum number of interpolation points (>n+1) (2*n+1)
     nfmax   [int] Maximum number of function evaluations (>n+1) (100)
     gtol    [dbl] Tolerance for the 2-norm of the model gradient (1e-4)
     delta   [dbl] Positive trust region radius (.1)
@@ -82,7 +82,7 @@ def pounders(fun, X0, n, mpmax, nfmax, gtol, delta, nfs, m, F0, xkin, L, U, prin
     # elif spsolver == 3:
     #     from minq8 import minq8
 
-    [flag, X0, mpmax, F0, L, U, xkin] = checkinputss(fun, X0, n, mpmax, nfmax, gtol, delta, nfs, m, F0, xkin, L, U)
+    [flag, X0, npmax, F0, L, U, xkin] = checkinputss(fun, X0, n, npmax, nfmax, gtol, delta, nfs, m, F0, xkin, L, U)
     if flag == -1:
         X = []
         F = []
@@ -132,7 +132,7 @@ def pounders(fun, X0, n, mpmax, nfmax, gtol, delta, nfs, m, F0, xkin, L, U, prin
         for i in range(nf + 1):
             D = X[i] - X[xkin]
             Res[i, :] = (F[i, :] - Cres) - 0.5 * D @ np.tensordot(D.T, Hres, 1)
-        [Mdir, mp, valid, Gres, Hresdel, Mind] = formquad(X[0 : nf + 1, :], Res[0 : nf + 1, :], delta, xkin, mpmax, Par, 0)
+        [Mdir, mp, valid, Gres, Hresdel, Mind] = formquad(X[0 : nf + 1, :], Res[0 : nf + 1, :], delta, xkin, npmax, Par, 0)
         if mp < n:
             [Mdir, mp] = bmpts(X[xkin], Mdir[0 : n - mp, :], L, U, delta, Par[2])
             for i in range(int(min(n - mp, nfmax - (nf + 1)))):
@@ -149,7 +149,7 @@ def pounders(fun, X0, n, mpmax, nfmax, gtol, delta, nfs, m, F0, xkin, L, U, prin
                 Res[nf, :] = (F[nf, :] - Cres) - 0.5 * D @ np.tensordot(D.T, Hres, 1)
             if nf + 1 >= nfmax:
                 break
-            [_, mp, valid, Gres, Hresdel, Mind] = formquad(X[0 : nf + 1, :], Res[0 : nf + 1, :], delta, xkin, mpmax, Par, False)
+            [_, mp, valid, Gres, Hresdel, Mind] = formquad(X[0 : nf + 1, :], Res[0 : nf + 1, :], delta, xkin, npmax, Par, False)
             if mp < n:
                 X, F, flag = prepare_outputs_before_return(X, F, nf, -5)
                 return
@@ -179,7 +179,7 @@ def pounders(fun, X0, n, mpmax, nfmax, gtol, delta, nfs, m, F0, xkin, L, U, prin
         # 2. Critically test invoked if the projected model gradient is small
         if ng < gtol:
             delta = max(gtol, np.max(np.abs(X[xkin])) * eps)
-            [Mdir, _, valid, _, _, _] = formquad(X[: nf + 1, :], F[: nf + 1, :], delta, xkin, mpmax, Par, 1)
+            [Mdir, _, valid, _, _, _] = formquad(X[: nf + 1, :], F[: nf + 1, :], delta, xkin, npmax, Par, 1)
             if not valid:
                 [Mdir, mp] = bmpts(X[xkin], Mdir, L, U, delta, Par[2])
                 for i in range(min(n - mp, nfmax - (nf + 1))):
@@ -195,7 +195,7 @@ def pounders(fun, X0, n, mpmax, nfmax, gtol, delta, nfs, m, F0, xkin, L, U, prin
                 if nf + 1 >= nfmax:
                     break
                 # Recalculate gradient based on a MFN model
-                [_, _, valid, Gres, Hres, Mind] = formquad(X[: nf + 1, :], F[: nf + 1, :], delta, xkin, mpmax, Par, 0)
+                [_, _, valid, Gres, Hres, Mind] = formquad(X[: nf + 1, :], F[: nf + 1, :], delta, xkin, npmax, Par, 0)
                 G, H = combinemodels(Cres, Gres, Hres)
                 ind_Lnotbinding = (X[xkin] > L) * (G.T > 0)
                 ind_Unotbinding = (X[xkin] < U) * (G.T < 0)
@@ -271,13 +271,13 @@ def pounders(fun, X0, n, mpmax, nfmax, gtol, delta, nfs, m, F0, xkin, L, U, prin
         # 5. Evaluate a model-improving point if necessary
         if not valid and (nf + 1 < nfmax) and (rho < eta1):  # Implies xkin, delta unchanged
             # Need to check because model may be valid after Xsp evaluation
-            [Mdir, mp, valid, _, _, _] = formquad(X[: nf + 1, :], F[: nf + 1, :], delta, xkin, mpmax, Par, 1)
+            [Mdir, mp, valid, _, _, _] = formquad(X[: nf + 1, :], F[: nf + 1, :], delta, xkin, npmax, Par, 1)
             if not valid:  # ! One strategy for choosing model-improving point:
                 # Update model (exists because delta & xkin unchanged)
                 for i in range(nf + 1):
                     D = X[i, :] - X[xkin]
                     Res[i, :] = (F[i, :] - Cres) - 0.5 * D @ np.tensordot(D.T, Hres, 1)
-                [_, _, valid, Gres, Hresdel, Mind] = formquad(X[: nf + 1, :], Res[: nf + 1, :], delta, xkin, mpmax, Par, False)
+                [_, _, valid, Gres, Hresdel, Mind] = formquad(X[: nf + 1, :], Res[: nf + 1, :], delta, xkin, npmax, Par, False)
                 Hres = Hres + Hresdel
                 # Update for modelimp; Cres unchanged b/c xkin unchanged
                 G, H = combinemodels(Cres, Gres, Hres)

pounders/py/tests/old_unit_tests/test_formquad.py

@@ -11,7 +11,7 @@ def test_formquad1(self):
         F = dictionaryData["F"]
         delta = dictionaryData["delta"]
         xkin = dictionaryData["xkin"]
-        mpmax = dictionaryData["mpmax"]
+        npmax = dictionaryData["npmax"]
         Pars = dictionaryData["Pars"]
         vf = dictionaryData["vf"]
         Mdir = dictionaryData["Mdir"]
@@ -22,7 +22,7 @@ def test_formquad1(self):
         G = dictionaryData["G"]
         H = dictionaryData["H"]
         Mind = dictionaryData["Mind"]
-        [MdirOut, mpOut, validOut, GOut, HOut, MindOut] = formquad(X, F, delta, xkin, mpmax, Pars, vf)
+        [MdirOut, mpOut, validOut, GOut, HOut, MindOut] = formquad(X, F, delta, xkin, npmax, Pars, vf)
         assert mpOut == mp
         assert np.shape(GOut) == np.shape(G)
         assert np.shape(HOut) == np.shape(H)
@@ -37,7 +37,7 @@ def test_formquad2(self):
         F = dictionaryData["F"]
         delta = dictionaryData["delta"]
         xkin = dictionaryData["xkin"]
-        mpmax = dictionaryData["mpmax"]
+        npmax = dictionaryData["npmax"]
         Pars = dictionaryData["Pars"]
         vf = dictionaryData["vf"]
         mp = dictionaryData["mp"]
@@ -47,7 +47,7 @@ def test_formquad2(self):
         Gres = dictionaryData["Gres"]
         Hresdel = dictionaryData["Hresdel"]
         Mind = dictionaryData["Mind"]
-        [_, mpOut, validOut, GOut, HOut, MindOut] = formquad(X, F, delta, xkin, mpmax, Pars, vf)
+        [_, mpOut, validOut, GOut, HOut, MindOut] = formquad(X, F, delta, xkin, npmax, Pars, vf)
         assert mpOut == mp
         assert np.linalg.norm(Gres - GOut) < 10**-10
         assert np.linalg.norm(Hresdel - HOut) < 10**-10
@@ -61,7 +61,7 @@ def test_formquad3(self):
         F = dictionaryData["F"]
         delta = dictionaryData["delta"]
         xkin = dictionaryData["xkin"]
-        mpmax = dictionaryData["mpmax"]
+        npmax = dictionaryData["npmax"]
         Pars = dictionaryData["Pars"]
         vf = dictionaryData["vf"]
         mp = dictionaryData["mp"]
@@ -72,7 +72,7 @@ def test_formquad3(self):
         H = dictionaryData["H"]
         Mind = dictionaryData["Mind"]
         Mdir = dictionaryData["Mdir"]
-        [MdirOut, mpOut, validOut, GOut, HOut, MindOut] = formquad(X, F, delta, xkin, mpmax, Pars, vf)
+        [MdirOut, mpOut, validOut, GOut, HOut, MindOut] = formquad(X, F, delta, xkin, npmax, Pars, vf)
         assert mpOut == mp
         assert np.linalg.norm(G - GOut) < 10**-10
         assert np.linalg.norm(H - HOut) < 10**-10
@@ -87,7 +87,7 @@ def test_formquad4(self):
         F = dictionaryData["F"]
         delta = dictionaryData["delta"]
         xkin = dictionaryData["xkin"]
-        mpmax = dictionaryData["mpmax"]
+        npmax = dictionaryData["npmax"]
         Pars = dictionaryData["Pars"]
         vf = dictionaryData["vf"]
         mp = dictionaryData["mp"]
@@ -98,7 +98,7 @@ def test_formquad4(self):
         H = dictionaryData["H"]
         Mind = dictionaryData["Mind"]
         Mdir = dictionaryData["Mdir"]
-        [MdirOut, mpOut, validOut, GOut, HOut, MindOut] = formquad(X, F, delta, xkin, mpmax, Pars, vf)
+        [MdirOut, mpOut, validOut, GOut, HOut, MindOut] = formquad(X, F, delta, xkin, npmax, Pars, vf)
         assert mpOut == mp - 1
         assert np.linalg.norm(G - GOut) < 10**-10
         assert np.linalg.norm(H - HOut) < 10**-10