Fix high-dimensional es-search (#38)

* fix: fixed High dimensions broadcasting operation issue

* add: added HD test for pybads

* chore: suppress warning test search

pybads/search/es_search.py

@@ -310,13 +310,9 @@ def ucov(U, u, w, ub, lb, scale, periodic_vars=None):
     u_shift = U_tmp - u_tmp
 
     if w.size != 0:
-        weights = w.reshape(
-            -1, *[1] * U.shape[1]
-        )  # For broadcasting weighted sum
-        C = np.sum(weights * (u_shift.T @ u_shift), axis=0)
-        C = C.reshape(
-            (U.shape[1], U.shape[1])
-        )  # Remove the extra dimesion from the broadcast result
+        weights = w.reshape(-1, *([1] * u_shift.ndim))  # For broadcasting weighted sum
+        C = np.matmul(u_shift.transpose(), weights * u_shift)
+        C = np.sum(C, axis=0)
     else:
         C = u_shift.T @ u_shift
 

pybads/testing/bads/scripts/bads_quadratic_noisy_example.py

@@ -3,11 +3,12 @@
 from pybads.bads import BADSDump
 from pybads.function_examples import quadratic_unknown_noisy_fcn, extra_noisy_quadratic_fcn, quadratic_hetsk_noisy_fcn, rosebrocks_hetsk_noisy_fcn
 
-x0 = np.array([[-3, -3]]);        # Starting point
-lb = np.array([[-5, -5]])     # Lower bounds
-ub = np.array([[5, 5]])       # Upper bounds
-plb = np.array([[-2, -2]])      # Plausible lower bounds
-pub = np.array([[2, 2]])        # Plausible upper bounds
+D = 3
+x0 = np.ones((1, D)) * -3      # Starting point
+lb = np.ones((1, D)) * -7      # Lower bounds
+ub = np.ones((1, D)) * 7       # Upper bounds
+plb = np.ones((1, D)) * -3     # Plausible lower bounds
+pub = np.ones((1, D)) * 3     # Plausible upper bounds
 
 title = 'Noise objective function'
 print("\n *** Example 3: " + title)

pybads/testing/bads/search/test_search.py

@@ -66,7 +66,7 @@ def test_search():
     pub = np.array([[5, 5, 5]])  # Plausible upper bounds
     D = 3
     bads = BADS(rosenbrocks_fcn, x0, lb, ub, plb, pub)
-    bads.options["fun_eval_start"] = 0
+    bads.options["fun_eval_start"] = 10
     gp, Ns_gp, sn2hpd, hyp_dict = bads._init_optimization_()
 
     es_iter = bads.options["n_search_iter"]
@@ -112,7 +112,7 @@ def test_search_hedge():
     D = 3
 
     bads = BADS(rosenbrocks_fcn, x0, lb, ub, plb, pub)
-    bads.options["fun_eval_start"] = 0
+    bads.options["fun_eval_start"] = 10
     gp, Ns_gp, sn2hpd, hyp_dict = bads._init_optimization_()
 
     search_hedge = ESSearchHedge(bads.options["search_method"], bads.options)
@@ -124,7 +124,6 @@ def test_search_hedge():
     assert us.size == 3 and (np.isscalar(z) or z.size == 1)
     assert np.all(gp.y >= z)
 
-
 def test_u_cov():
     U = np.array(
         [
@@ -142,7 +141,6 @@ def test_u_cov():
     C = ucov(U, u0, w, ub, lb, 1)
     assert C.shape == (U.shape[1], U.shape[1])
 
-
 def test_grid_search_neighbors():
     x0 = np.array([[0, 0]])
     # Starting point
@@ -153,7 +151,7 @@ def test_grid_search_neighbors():
     D = 2
 
     bads = BADS(rosenbrocks_fcn, x0, lb, ub, plb, pub)
-    bads.options["fun_eval_start"] = 0
+    bads.options["fun_eval_start"] = 10
     gp, Ns_gp, sn2hpd, hyp_dict = bads._init_optimization_()
     gp.X = np.array([[0, 0], [-0.1055, 0.4570], [-0.3555, -0.7930]])
     f = FunctionLogger(rosenbrocks_fcn, D, False, 0)
@@ -175,3 +173,4 @@ def test_grid_search_neighbors():
         and np.isclose(result[1, 0], -0.1055, 1e-3)
         and np.isclose(result[2, 0], -0.3555, 1e-3)
     )
+

pybads/testing/bads/test_bads_optimization.py

@@ -6,8 +6,7 @@
 from pybads.bads import BADS
 
 
-def get_test_opt_conf():
-    D = 3
+def get_test_opt_conf(D=3):
     x0 = np.ones((1, D)) * 4
     LB = -100*np.ones(D)                # Lower bound
     UB = 100*np.ones(D)                 # Upper bound
@@ -17,18 +16,19 @@ def get_test_opt_conf():
     return D, x0, LB, UB, PLB, PUB, tol_errs
 
 def run_bads(fun, x0, LB, UB, PLB, PUB, tol_errs, f_min,
-            oracle_fun=None, non_box_cons=None, he_noise_flag=False,
-            assert_flag=False):
+            oracle_fun=None, non_box_cons=None, uncertainty_handling=False,
+            assert_flag=False, max_fun_evals=None):
     options = {}
     options["display"] = "full"#debug_flag = True
 
-    if he_noise_flag:
-        options["specify_target_noise"] = True
+    if uncertainty_handling > 0:
         options["uncertainty_handling"] = True
-        #options["noise_final_samples"] = 100
-        options["max_fun_evals"] = 200
+        options["max_fun_evals"] = 200 if max_fun_evals is None  else max_fun_evals
+        options["specify_target_noise"] = True
+        if uncertainty_handling > 1:
+            options["specify_target_noise"] = True
     else:
-        options["max_fun_evals"] = 100
+        options["max_fun_evals"] = 100 if max_fun_evals is None else max_fun_evals
 
     optimize_result = BADS(fun=fun, x0=x0, lower_bounds=LB,
                                     upper_bounds=UB, plausible_lower_bounds=PLB,
@@ -47,11 +47,17 @@ def run_bads(fun, x0, LB, UB, PLB, PUB, tol_errs, f_min,
     if assert_flag:
         assert np.any(err < tol_errs), f"Error {err} is not smaller than tolerance {tol_errs} when optimizing {fun.__name__}."
 
+    return optimize_result, err
+
 def test_ellipsoid_opt():
     D, x0, LB, UB, PLB, PUB, tol_errs = get_test_opt_conf()
     fun = lambda x: np.sum((np.atleast_2d(x) / np.arange(1, len(x) + 1) ** 2) ** 2)
     run_bads(fun, x0, LB, UB, PLB, PUB, tol_errs, f_min=0.0, assert_flag=True)
-
+
+def test_high_dim_opt():
+    D, x0, LB, UB, PLB, PUB, tol_errs = get_test_opt_conf(D=60)
+    fun = lambda x: np.sum((np.atleast_2d(x) / np.arange(1, len(x) + 1) ** 2) ** 2)
+    run_bads(fun, x0, LB, UB, PLB, PUB, tol_errs, f_min=0.0, assert_flag=False, max_fun_evals=200)
 
 def test_sphere_opt():
     D, x0, LB, UB, PLB, PUB, tol_errs = get_test_opt_conf()
@@ -77,4 +83,4 @@ def test_he_noisy_sphere_opt():
     D, x0, LB, UB, PLB, PUB, tol_errs = get_test_opt_conf()
     fun = he_noisy_sphere
     oracle_fun = lambda x: np.sum(np.atleast_2d(x)**2, axis=1)                          # True objective function
-    run_bads(fun, x0, LB, UB, PLB, PUB, tol_errs, f_min=0.0, oracle_fun=oracle_fun, he_noise_flag=True, assert_flag=True)
+    run_bads(fun, x0, LB, UB, PLB, PUB, tol_errs, f_min=0.0, oracle_fun=oracle_fun, uncertainty_handling=2, assert_flag=True)