ocba var handling

bartzbeielstein · bartzbeielstein · commit bad219f51909 · 2024-01-01T10:09:23.000+01:00
diff --git a/notebooks/testKriging.ipynb b/notebooks/testKriging.ipynb
@@ -1392,12 +1392,161 @@
     "S.min_y"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## test ocba"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {},
    "outputs": [],
-   "source": []
+   "source": [
+    "\n",
+    "import copy\n",
+    "import numpy as np\n",
+    "from spotPython.fun.objectivefunctions import analytical\n",
+    "from spotPython.spot import spot\n",
+    "from spotPython.budget.ocba import get_ocba\n",
+    "\n",
+    "# Test based on the example from the book:\n",
+    "# Chun-Hung Chen and Loo Hay Lee:\n",
+    "#     Stochastic Simulation Optimization: An Optimal Computer Budget Allocation,\n",
+    "#     pp. 49 and pp. 215\n",
+    "#     p. 49:\n",
+    "#     mean_y = np.array([1,2,3,4,5])\n",
+    "#     var_y = np.array([1,1,9,9,4])\n",
+    "#     get_ocba(mean_y, var_y, 50)\n",
+    "#     [11  9 19  9  2]\n",
+    "\n",
+    "fun = analytical().fun_linear\n",
+    "fun_control = {\"sigma\": 0.001,\n",
+    "               \"seed\": 123}\n",
+    "spot_1_noisy = spot.Spot(fun=fun,\n",
+    "                   lower = np.array([-1]),\n",
+    "                   upper = np.array([1]),\n",
+    "                   fun_evals = 20,\n",
+    "                   fun_repeats = 2,\n",
+    "                   noise = True,\n",
+    "                   ocba_delta=1,\n",
+    "                   seed=123,\n",
+    "                   show_models=False,\n",
+    "                   fun_control = fun_control,\n",
+    "                   design_control={\"init_size\": 3,\n",
+    "                                   \"repeats\": 2},\n",
+    "                   surrogate_control={\"noise\": True})\n",
+    "spot_1_noisy.run()\n",
+    "spot_2 = copy.deepcopy(spot_1_noisy)\n",
+    "spot_2.mean_y = np.array([1,2,3,4,5])\n",
+    "spot_2.var_y = np.array([1,1,9,9,4])\n",
+    "n = 50\n",
+    "o = get_ocba(spot_2.mean_y, spot_2.var_y, n)\n",
+    "assert sum(o) == 50\n",
+    "assert (o == np.array([[11, 9, 19, 9, 2]])).all()\n",
+    "o"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "Seed set to 123\n"
+     ]
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "S.X: [[ 0.          1.        ]\n",
+      " [ 1.          0.        ]\n",
+      " [ 1.          1.        ]\n",
+      " [ 1.          1.        ]\n",
+      " [ 0.54509876 -0.36921401]\n",
+      " [ 0.54509876 -0.36921401]\n",
+      " [ 0.18642675  0.87708546]\n",
+      " [ 0.18642675  0.87708546]\n",
+      " [-0.45060393 -0.208063  ]\n",
+      " [-0.45060393 -0.208063  ]]\n",
+      "S.y: [0.98021757 0.99264427 2.02575851 2.00387949 0.45185626 0.44499372\n",
+      " 0.79130456 0.81487288 0.24000221 0.23988634]\n",
+      "X_shape_before: (10, 2)\n",
+      "y_size_before: 10\n",
+      "S.X: [[ 0.          1.        ]\n",
+      " [ 1.          0.        ]\n",
+      " [ 1.          1.        ]\n",
+      " [ 1.          1.        ]\n",
+      " [ 0.54509876 -0.36921401]\n",
+      " [ 0.54509876 -0.36921401]\n",
+      " [ 0.18642675  0.87708546]\n",
+      " [ 0.18642675  0.87708546]\n",
+      " [-0.45060393 -0.208063  ]\n",
+      " [-0.45060393 -0.208063  ]\n",
+      " [-0.39841465 -0.21105872]\n",
+      " [-0.39841465 -0.21105872]]\n",
+      "S.y: [0.98021757 0.99264427 2.02575851 2.00387949 0.45185626 0.44499372\n",
+      " 0.79130456 0.81487288 0.24000221 0.23988634 0.18349759 0.19592429]\n",
+      "S.n_points: 1\n",
+      "S.ocba_delta: 1\n",
+      "X_shape_after: (12, 2)\n",
+      "y_size_after: 12\n"
+     ]
+    }
+   ],
+   "source": [
+    "import numpy as np\n",
+    "from spotPython.fun.objectivefunctions import analytical\n",
+    "from spotPython.spot import spot\n",
+    "from spotPython.utils.init import fun_control_init\n",
+    "# number of initial points:\n",
+    "ni = 3\n",
+    "X_start = np.array([[0, 1], [1, 0], [1, 1], [1, 1]])\n",
+    "\n",
+    "fun = analytical().fun_sphere\n",
+    "fun_control = fun_control_init(\n",
+    "        sigma=0.02,\n",
+    "        seed=123,)\n",
+    "lower = np.array([-1, -1])\n",
+    "upper = np.array([1, 1])\n",
+    "design_control={\"init_size\": ni,\n",
+    "                \"repeats\": 2}\n",
+    "\n",
+    "S = spot.Spot(fun=fun,\n",
+    "              noise=True,\n",
+    "              fun_repeats=2,\n",
+    "              n_points=1,\n",
+    "              ocba_delta=1,\n",
+    "               log_level = 10,\n",
+    "            lower = lower,\n",
+    "            upper= upper,\n",
+    "            show_progress=False,\n",
+    "            design_control=design_control,\n",
+    "            fun_control=fun_control\n",
+    ")\n",
+    "S.initialize_design(X_start=X_start)\n",
+    "print(f\"S.X: {S.X}\")\n",
+    "print(f\"S.y: {S.y}\")\n",
+    "X_shape_before = S.X.shape\n",
+    "print(f\"X_shape_before: {X_shape_before}\")\n",
+    "print(f\"y_size_before: {S.y.size}\")\n",
+    "y_size_before = S.y.size\n",
+    "S.update_stats()\n",
+    "S.fit_surrogate()\n",
+    "S.update_design()\n",
+    "print(f\"S.X: {S.X}\")\n",
+    "print(f\"S.y: {S.y}\")\n",
+    "print(f\"S.n_points: {S.n_points}\")\n",
+    "print(f\"S.ocba_delta: {S.ocba_delta}\")\n",
+    "print(f\"X_shape_after: {S.X.shape}\")\n",
+    "print(f\"y_size_after: {S.y.size}\")\n"
+   ]
   },
   {
    "cell_type": "code",
diff --git a/pyproject.toml b/pyproject.toml
@@ -7,7 +7,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "spotPython"
-version = "0.8.15"
+version = "0.8.16"
 authors = [
   { name="T. Bartz-Beielstein", email="tbb@bartzundbartz.de" }
 ]
diff --git a/src/spotPython/budget/ocba.py b/src/spotPython/budget/ocba.py
@@ -5,10 +5,10 @@
 from spotPython.utils.aggregate import get_ranks
 from numpy import int32, float64
 from numpy import argpartition, repeat
-from numpy import zeros, square, sqrt, full, around
+from numpy import zeros, square, sqrt, full, around, array
 
 
-def get_ocba(means, vars, delta) -> int32:
+def get_ocba(means, vars, delta, verbose=False) -> array:
     """
     Optimal Computer Budget Allocation (OCBA)
 
@@ -25,9 +25,14 @@ def get_ocba(means, vars, delta) -> int32:
         https://github.com/TomMonks/sim-tools
 
     Args:
-        means (numpy.array): An array of means.
-        vars (numpy.array): An array of variances.
-        delta (int): The incremental budget.
+        means (numpy.array):
+            An array of means.
+        vars (numpy.array):
+            An array of variances.
+        delta (int):
+            The incremental budget.
+        verbose (bool):
+            If True, print the results.
 
     Returns:
         (numpy.array): An array of budget recommendations.
@@ -36,11 +41,11 @@ def get_ocba(means, vars, delta) -> int32:
         The implementation is based on the pseudo-code in the Chen et al. (p. 49), see [1].
 
     Examples:
-        From the Chen et al. book (p. 49):
+        >>> # From the Chen et al. book (p. 49):
         >>> mean_y = np.array([1,2,3,4,5])
             var_y = np.array([1,1,9,9,4])
             get_ocba(mean_y, var_y, 50)
-            [11  9 19  9  2]
+            array([11,  9, 19,  9,  2])
     """
     n_designs = means.shape[0]
     allocations = zeros(n_designs, int32)
@@ -58,6 +63,23 @@ def get_ocba(means, vars, delta) -> int32:
     more_runs = full(n_designs, True, dtype=bool)
     add_budget = zeros(n_designs, dtype=float)
     more_alloc = True
+    if verbose:
+        print("\nIn get_ocba():")
+        print(f"means: {means}")
+        print(f"vars: {vars}")
+        print(f"delta: {delta}")
+        print(f"n_designs: {n_designs}")
+        print(f"Allocations: {allocations}")
+        print(f"Ratios: {ratios}")
+        print(f"Budget: {budget}")
+        print(f"Ranks: {ranks}")
+        print(f"Best: {best}")
+        print(f"Second best: {second_best}")
+        print(f"Select: {select}")
+        print(f"Temp: {temp}")
+        print(f"More runs: {more_runs}")
+        print(f"Add budget: {add_budget}")
+        print(f"More allocations: {more_alloc}")
     while more_alloc:
         more_alloc = False
         ratio_s = (more_runs * ratios).sum()
@@ -66,6 +88,11 @@ def get_ocba(means, vars, delta) -> int32:
         mask = add_budget < allocations
         add_budget[mask] = allocations[mask]
         more_runs[mask] = 0
+        if verbose:
+            print("\nIn more_alloc:")
+            print(f"ratio_s: {ratio_s}")
+            print(f"more_runs: {more_runs}")
+            print(f"add_budget: {add_budget}")
         if mask.sum() > 0:
             more_alloc = True
         if more_alloc:
diff --git a/src/spotPython/spot/spot.py b/src/spotPython/spot/spot.py
@@ -249,7 +249,7 @@ def __init__(
         # specified once:
         if len(self.var_type) < self.k:
             self.var_type = self.var_type * self.k
-            logger.warning("Warning: All variable types forced to 'num'.")
+            logger.warning("All variable types forced to 'num'.")
         self.infill_criterion = infill_criterion
         # Bounds
         de_bounds = []
@@ -749,14 +749,19 @@ def update_design(self) -> None:
         """
         # OCBA (only if noise). Determination of the OCBA points depends on the
         # old X and y values.
-        if self.noise and self.ocba_delta > 0:
+        if self.noise and self.ocba_delta > 0 and not np.all(self.var_y > 0):
+            logger.warning(
+                "self.var_y <= 0. OCBA points are not generated, because there are points with no variance information."
+            )
+            logger.debug("In update_design(): self.var_y: %s", self.var_y)
+        if self.noise and self.ocba_delta > 0 and np.all(self.var_y > 0):
             X_ocba = get_ocba_X(self.mean_X, self.mean_y, self.var_y, self.ocba_delta)
         else:
             X_ocba = None
         # Determine the new X0 values based on the old X and y values:
         X0 = self.get_new_X0()
         # Append OCBA points to the new design points:
-        if self.noise and self.ocba_delta > 0:
+        if self.noise and self.ocba_delta > 0 and np.all(self.var_y > 0):
             X0 = append(X_ocba, X0, axis=0)
         X_all = self.to_all_dim_if_needed(X0)
         logger.debug(
diff --git a/test/test_update_design.py b/test/test_update_design.py
@@ -31,7 +31,7 @@ def test_update_design():
 def test_update_design_with_repeats_and_ocba():
     # number of initial points:
     ni = 3
-    X_start = np.array([[0, 1], [1, 0], [1, 1], [1, 1]])
+    X_start = np.array([[1, 0], [1, 0], [1, 1], [1, 1]])
     fun = analytical().fun_sphere
     fun_control = fun_control_init(
             sigma=0.02,
@@ -61,3 +61,39 @@ def test_update_design_with_repeats_and_ocba():
     assert X_shape_before[0] + S.n_points * S.fun_repeats + S.ocba_delta == S.X.shape[0]
     assert X_shape_before[1] == S.X.shape[1]
     assert y_size_before + S.n_points * S.fun_repeats + S.ocba_delta == S.y.size
+
+def test_update_design_with_repeats_and_ocba_no_var():
+    # number of initial points:
+    ni = 3
+    # The first two X_start points have no repeats and therefore no variance
+    X_start = np.array([[0, 1], [1, 0], [1, 1], [1, 1]])
+    fun = analytical().fun_sphere
+    fun_control = fun_control_init(
+            sigma=0.02,
+            seed=123,)
+    lower = np.array([-1, -1])
+    upper = np.array([1, 1])
+    design_control={"init_size": ni,
+                    "repeats": 2}
+    S = spot.Spot(fun=fun,
+                noise=True,
+                fun_repeats=2,
+                n_points=1,
+                ocba_delta=1,
+                lower = lower,
+                upper= upper,
+                show_progress=True,
+                design_control=design_control,
+                fun_control=fun_control
+    )
+    S.initialize_design(X_start=X_start)
+    X_shape_before = S.X.shape
+    y_size_before = S.y.size
+    S.update_stats()
+    S.fit_surrogate()
+    S.update_design()
+    # compare the shapes of the X and y values before and after the update_design method
+    # no ocba points are added because the first two points have no variance
+    assert X_shape_before[0] + S.n_points * S.fun_repeats == S.X.shape[0]
+    assert X_shape_before[1] == S.X.shape[1]
+    assert y_size_before + S.n_points * S.fun_repeats == S.y.size

Original file line number	Diff line number	Diff line change
`@@ -7,7 +7,7 @@ build-backend = "setuptools.build_meta"`
`7`	`7`
`8`	`8`	`[project]`
`9`	`9`	`name = "spotPython"`
`10`		`-version = "0.8.15"`
	`10`	`+version = "0.8.16"`
`11`	`11`	`authors = [`
`12`	`12`	`{ name="T. Bartz-Beielstein", email="tbb@bartzundbartz.de" }`
`13`	`13`	`]`