New/py goombah

goombah/py/goombah.py

@@ -0,0 +1,172 @@
+# This code solves the problem
+#     minimize h(F(x))
+# where x is an [n by 1] vector, F is a blackbox function mapping from R^n to
+# R^p, and h is a nonsmooth function mapping from R^p to R.
+#
+#
+# Inputs:
+#  hfun:    [func handle] Evaluates h, returning the [scalar] function
+#                         value and [k x m] subgradients for all k limiting
+#                         gradients at the point given.
+#  Ffun:    [func handle] Evaluates F, the black box simulation, returning
+#                         a [1 x m] vector.
+#  nfmax:   [int]         Maximum number of function evaluations.
+#  x0:      [1 x n dbl]   Starting point.
+#  L:       [1 x n dbl]   Lower bounds.
+#  U:       [1 x n dbl]   Upper bounds.
+#  GAMS_options:
+#  subprob_switch:
+#
+# Outputs:
+#   X:      [nfmax x n]   Points evaluated
+#   F:      [nfmax x p]   Their simulation values
+#   h:      [nfmax x 1]   The values h(F(x))
+#   xkin:   [int]         Current trust region center
+
+import numpy as np
+import sys
+from check_inputs_and_initialize import check_inputs_and_initialize
+from call_user_scripts import call_user_scripts
+from checkinputss import checkinputss
+from build_p_models import build_p_models
+from save_quadratics_call_GAMS import save_quadratics_call_GAMS
+from choose_generator_set import choose_generator_set
+from minimize_affine_envelope import minimize_affine_envelope
+
+
+def goombah(hfun, Ffun, nfmax, x0, L, U, GAMS_options, subprob_switch):
+    try:
+        F0 = Ffun(x0)
+        F0 = F0.flatten()
+    except:
+        print("Problem using Ffun. Exiting")
+        X, F, h, xkin = [], [], [], []
+        flag = -1
+        return X, F, h, xkin, flag
+
+    n, delta, _, fq_pars, tol, X, F, h, Hash, nf, successful, xkin, Hres = check_inputs_and_initialize(x0, F0, nfmax)
+    flag, x0, __, F0, L, U, xkin = checkinputss(hfun, x0, n, fq_pars["npmax"], nfmax, tol["gtol"], delta, 1, len(F0), F0, xkin, L, U)
+
+    h[nf], _, hashes_at_nf = hfun(F[nf, :])
+    Hash[nf] = hashes_at_nf
+
+    I_n = np.eye(n)
+    for i in range(n):
+        nf, X, F, h, Hash, _ = call_user_scripts(nf, X, F, h, Hash, Ffun, hfun, X[xkin] + delta * I_n[i, :], tol, L, U, 1)
+
+    H_mm = np.zeros((n, n))
+    beta_exp = 1.0
+
+    while nf + 1 < nfmax and delta > tol["mindelta"]:
+        xkin = int(np.argmin(h[: nf + 1]))
+
+        Gres, Hres, X, F, h, nf, Hash = build_p_models(nf, nfmax, xkin, delta, F, X, h, Hres, fq_pars, tol, hfun, Ffun, Hash, L, U)
+
+        if len(Gres) == 0:
+            print("Empty Gres. Delta =", delta)
+            X = X[: nf + 1]
+            F = F[: nf + 1]
+            h = h[: nf + 1]
+            flag = -1
+            return X, F, h, xkin, flag
+
+        if nf + 1 >= nfmax:
+            flag = 0
+            return X, F, h, xkin, flag
+
+        Low = np.maximum(L - X[xkin], -delta)
+        Upp = np.minimum(U - X[xkin], delta)
+
+        sk, pred_dec = save_quadratics_call_GAMS(Hres, Gres, F[xkin, :], Low, Upp, X[xkin], X[xkin], h[xkin], GAMS_options, hfun)
+
+        if pred_dec == 0:
+            rho_k = -np.inf
+        else:
+            nf, X, F, h, Hash, _ = call_user_scripts(nf, X, F, h, Hash, Ffun, hfun, X[xkin] + sk, tol, L, U, 1)
+            rho_k = (h[xkin] - h[nf]) / min(1.0, delta ** (1.0 + beta_exp))
+
+        if rho_k > tol["eta1"]:
+            if np.linalg.norm(X[xkin] - X[nf], np.inf) >= 0.8 * delta:
+                delta = delta * tol["gamma_inc"]
+            xkin = nf
+        else:
+            # Need to do one MS-P loop
+
+            bar_delta = delta
+
+            while nf + 1 < nfmax:
+                Gres, Hres, X, F, h, nf, Hash = build_p_models(nf, nfmax, xkin, delta, F, X, h, Hres, fq_pars, tol, hfun, Ffun, Hash, L, U)
+
+                if len(Gres) == 0:
+                    print("Model building failed. Empty Gres. Delta =", delta)
+                    X = X[: nf + 1]
+                    F = F[: nf + 1]
+                    h = h[: nf + 1]
+                    flag = -1
+                    return X, F, h, xkin, flag
+
+                if nf + 1 >= nfmax:
+                    flag = 0
+                    return X, F, h, xkin, flag
+
+                D_k, Act_Z_k, f_bar = choose_generator_set(X, Hash, 3, xkin, nf, delta, F, hfun)
+                G_k = np.dot(Gres, D_k)
+                beta = max(0, f_bar - h[xkin])
+
+                H_k = np.zeros((G_k.shape[1], n + 1, n + 1))
+                for i in range(G_k.shape[1]):
+                    for j in range(Hres.shape[2]):
+                        H_k[i, 1:, 1:] += D_k[j, i] * Hres[:, :, j]
+
+                Low = np.maximum(L - X[xkin], -delta)
+                Upp = np.minimum(U - X[xkin], delta)
+
+                s_k, tau_k, _, _ = minimize_affine_envelope(h[xkin], f_bar, beta, G_k, H_mm, delta, Low, Upp, H_k, subprob_switch)
+
+                Low = np.maximum(L - X[xkin], -1.0)
+                Upp = np.minimum(U - X[xkin], 1.0)
+                _, _, chi_k, _ = minimize_affine_envelope(h[xkin], f_bar, beta, G_k, np.zeros((n, n)), delta, Low, Upp, np.zeros((G_k.shape[1], n + 1, n + 1)), subprob_switch)
+
+                if chi_k <= tol["gtol"] and delta <= tol["mindelta"]:
+                    print("Convergence satisfied: small stationary measure and small delta")
+                    X = X[: nf + 1]
+                    F = F[: nf + 1]
+                    h = h[: nf + 1]
+                    flag = chi_k
+                    return X, F, h, xkin, flag
+
+                nf, X, F, h, Hash, hashes_at_nf = call_user_scripts(nf, X, F, h, Hash, Ffun, hfun, X[xkin] + s_k.T, tol, L, U, 1)
+
+                ared = h[xkin] - h[nf]
+                pred = -tau_k
+                rho_k = ared / pred
+
+                if rho_k >= tol["eta1"] and pred > 0:
+                    successful = True
+                    break
+                else:
+                    tmp_Act_Z_k = choose_generator_set(X, Hash, 3, xkin, nf, delta, F, hfun)[1]
+
+                    if all(item in tmp_Act_Z_k for item in Act_Z_k):
+                        if any(item in hashes_at_nf for item in Act_Z_k):
+                            successful = False
+                            break
+                        else:
+                            delta = tol["gamma_dec"] * delta
+            if successful:
+                xkin = nf
+                if rho_k > tol["eta3"] and np.linalg.norm(s_k, np.inf) > 0.8 * bar_delta:
+                    delta = bar_delta * tol["gamma_inc"]
+            else:
+                delta = max(bar_delta * tol["gamma_dec"], tol["mindelta"])
+
+        print(f"nf: {nf:8d}; fval: {h[xkin,0]:8e}; radius: {delta:8e};")
+
+    if nf + 1 >= nfmax:
+        flag = 0
+    else:
+        X = X[: nf + 1]
+        F = F[: nf + 1]
+        h = h[: nf + 1]
+
+    return X, F, h, xkin, flag

goombah/py/save_quadratics_call_GAMS.py

@@ -0,0 +1,65 @@
+# A general Python-to-GAMS interface.
+import numpy as np
+import subprocess
+
+
+def save_quadratics_call_GAMS(H, g, b, Low, Upp, x0, x1, val_at_x0, GAMS_options, hfun):
+    n, P = g.shape
+
+    # Loop over solvers
+    for i in GAMS_options["solvers"]:
+        solver_val = i
+
+        # Put problem data to a gdx file
+        # wgdx('quad_model_data', Ns, Ps, Hs_mod, gs_mod, bs_mod, x0s, x1s, solver, Lows, Upps)
+        # print('Data written to GDX file quads_model_data.gdx')
+
+        # # Copy the template gams file
+        # subprocess.run(['cp', GAMS_options['file'], f'./TRSP_{i}.gms'])
+
+        # # Perform the gams run
+        # flag = subprocess.call(['gams', f'TRSP_{i}.gms', 'lo=2'])
+
+        # assert flag == 0, f'gams run failed: rc = {flag}'
+        # assert os.path.exists(solGDX), f'Results file {solGDX} does not exist after gams run'
+
+        # print(f'TRSP_{i}.gms finished')
+
+        # # Now get the outputs from the GDX file produced by the GAMS run
+        # rs = {'name': 'modelStat', 'form': 'full'}
+        # r = rgdx(solGDX, rs)
+        # modelStat[i] = r['val']
+
+        # rs['name'] = 'solveStat'
+        # r = rgdx(solGDX, rs)
+        # solveStat[i] = r['val']
+
+        # rs['name'] = 'tau'
+        # rs['field'] = 'l'
+        # r = rgdx(solGDX, rs)
+        # obj_vals_GAMS[i] = r['val']
+
+        # rs['name'] = 'x'
+        # rs['uels'] = Ns['uels']
+        # r = rgdx(solGDX, rs)
+        # x = r['val']
+        # allx[i, :] = x
+
+    # Just randomly try values
+    num_rand_samp = 1000
+    obj_vals_PYTHON = np.zeros(num_rand_samp)
+    allx = np.random.uniform(x0 + Low, x0 + Upp, (num_rand_samp, n))
+    z = np.zeros(P)
+    for j, x in enumerate(allx):
+        z = 0.5 * (x - x0) @ ((x - x0) @ H) + (x - x0) @ g + b
+        obj_vals_PYTHON[j] = hfun(z)[0]
+
+    val_at_new = np.min(obj_vals_PYTHON)
+    ind = np.argmin(obj_vals_PYTHON)
+    s_k = allx[ind, :] - x0
+    pred_dec = val_at_x0 - val_at_new
+
+    if pred_dec < 0:
+        pred_dec = 0
+
+    return s_k, pred_dec

goombah/py/tests/compare_goombah_and_pounders.py

@@ -0,0 +1,85 @@
+import os
+import sys
+
+import numpy as np
+import scipy as sp
+
+sys.path.append("../../../../minq/py/minq5/")  # Needed for spsolver=2
+sys.path.append("../")
+from ibcdfo.pounders import pounders
+
+BenDFO_root = "../../../../BenDFO/"
+sys.path.append(BenDFO_root + "py/")  # Needed for spsolver=2
+sys.path.append("../../../manifold_sampling/py/")
+sys.path.append("../../../pounders/py/")
+sys.path.append("../../../manifold_sampling/py/h_examples")
+
+from sum_squared import sum_squared
+from dfoxs import dfoxs
+from calfun import calfun
+from goombah import goombah
+
+os.makedirs("benchmark_results", exist_ok=True)
+
+dfo = np.loadtxt(BenDFO_root + "data/dfo.dat")
+
+spsolver = 2  # TRSP solver
+nfmax = 50
+gtol = 1e-13
+factor = 10
+
+Results = {}
+GAMS_options = {}
+hfun = lambda F: np.sum(F**2)
+
+for row, (nprob, n, m, factor_power) in enumerate(dfo):
+    n = int(n)
+    m = int(m)
+
+    def objective(y):
+        # It is possible to have python use the same objective values via
+        # octave. This can be slow on some systems. To (for example)
+        # test difference between matlab and python, used the following
+        # line and add "from oct2py import octave" on a system with octave
+        # installed.
+        # out = octave.feval("calfun_wrapper", y, m, nprob, "smooth", [], 1, 1)
+        out = calfun(y, m, int(nprob), "smooth", 0, vecout=True)
+        assert len(out) == m, "Incorrect output dimension"
+        return np.squeeze(out)
+
+    X0 = dfoxs(n, nprob, int(factor**factor_power))
+    npmax = 2 * n + 1  # Maximum number of interpolation points [2*n+1]
+    L = -np.inf * np.ones((1, n))  # 1-by-n Vector of lower bounds [zeros(1, n)]
+    U = np.inf * np.ones((1, n))  # 1-by-n Vector of upper bounds [ones(1, n)]
+    nfs = 1
+    F0 = np.zeros((1, m))
+    F0[0] = objective(X0)
+    xind = 0
+    delta = 0.1
+    printf = False
+
+    filename = "./benchmark_results/poundersM_and_GOOMBAH_nfmax=" + str(nfmax) + "_gtol=" + str(gtol) + "_prob=" + str(row) + "_spsolver=" + str(spsolver) + ".mat"
+
+    Results[row] = {}
+    for method in [1, 2]:
+        if method == 1:
+            [X, F, flag, xk_best] = pounders(objective, X0, n, npmax, nfmax, gtol, delta, nfs, m, F0, xind, L, U, printf, spsolver)
+            Results[row]["alg"] = "pounders4py"
+        elif method == 2:
+            GAMS_options["solvers"] = range(4)
+            [X, F, h, xk_best, flag] = goombah(sum_squared, objective, nfmax, X0, L, U, GAMS_options, "linprog")
+            Results[row]["alg"] = "goombah"
+
+        evals = F.shape[0]
+        h = np.zeros(evals)
+
+        for i in range(evals):
+            h[i] = hfun(F[i, :])
+
+        Results[row]["problem"] = "problem " + str(row) + " from More/Wild"
+        Results[row]["Fvec"] = F
+        Results[row]["H"] = h
+        Results[row]["X"] = X
+
+
+np.save("Results.npy", Results)

manifold_sampling/py/choose_generator_set.py

@@ -12,16 +12,19 @@ def choose_generator_set(X, Hash, gentype, xkin, nf, delta, F, hfun):
 
     elif gentype == 3:
         hxkin, _ = hfun(F[xkin, :], Act_Z_k)
+        hxkin = np.atleast_1d(hxkin)
         for i in [ind for ind in range(nf) if ind != xkin]:
             Act_tmp = Hash[i]
             h_i, _ = hfun(F[xkin], Act_tmp)
+            h_i = np.atleast_1d(h_i)
             if np.linalg.norm(X[xkin] - X[i], ord=np.inf) <= delta * (1 + 1e-8) and h_i[0] <= hxkin[0]:
                 Act_Z_k = np.concatenate((Act_Z_k, Act_tmp))
             elif np.linalg.norm(X[xkin] - X[i], ord=np.inf) <= delta**2 * (1 + 1e-8) and h_i[0] > hxkin[0]:
                 Act_Z_k = np.concatenate((Act_Z_k, Act_tmp))
 
     Act_Z_k = np.unique(Act_Z_k)
     f_k, D_k = hfun(F[xkin], Act_Z_k)
+    f_k = np.atleast_1d(f_k)
     unique_indices = np.unique(D_k, axis=1, return_index=True)[1]
     D_k = D_k[:, unique_indices]
     Act_Z_k = Act_Z_k[unique_indices]