Merge pull request #23 from IOHprofiler/fixing-weights

Fixing weights

modcma/modularcmaes.py

@@ -82,7 +82,8 @@ def mutate(self) -> None:
 
         y = np.dot(self.parameters.B, self.parameters.D * z)
         x = self.parameters.m + (s * y)
-        x, n_out_of_bounds = correct_bounds(x, 
+        x, n_out_of_bounds = correct_bounds(
+            x, 
             self.parameters.ub,
             self.parameters.lb,
             self.parameters.bound_correction
@@ -375,7 +376,7 @@ def correct_bounds(
     """
     out_of_bounds = np.logical_or(x > ub, x < lb)
     n_out_of_bounds = out_of_bounds.max(axis=0).sum()
-    if n_out_of_bounds == 0 or correction_method == None:
+    if n_out_of_bounds == 0 or correction_method is None:
         return x, n_out_of_bounds
 
     try:
@@ -464,7 +465,7 @@ def evaluate_bbob(
     fitness_func = IOH_function(
         fid, dim, instance, target_precision=target_precision, suite="BBOB"
     )
-    
+
     if logging:
         data_location = data_folder if os.path.isdir(data_folder) else os.getcwd()
         logger = IOH_logger(data_location, f"{label}F{fid}_{dim}D")

modcma/parameters.py

@@ -25,6 +25,8 @@ class Parameters(AnnotatedStruct):
     ----------
     d: int
         The dimensionality of the problem
+    x0: np.ndarray
+        Initial guess of the population center of mass.
     target: float = -float("inf")
         The absolute target of the optimization problem
     budget: int = None
@@ -37,7 +39,7 @@ class Parameters(AnnotatedStruct):
         The number of offspring in the population
     mu: int = None
         The number of parents in the population
-    init_sigma: float = .5
+    sigma0: float = .5
         The initial value of sigma (step size)
     a_tpa: float = .5
         Parameter used in TPA
@@ -236,12 +238,13 @@ class Parameters(AnnotatedStruct):
     """
 
     d: int
+    x0: np.ndarray = None
     target: float = -float("inf")
     budget: int = None
     n_generations: int = None
     lambda_: int = None
     mu: int = None
-    init_sigma: float = 0.5
+    sigma0: float = 0.5
     a_tpa: float = 0.5
     b_tpa: float = 0.0
     cs: float = None
@@ -346,6 +349,10 @@ def init_fixed_parameters(self) -> None:
         self.bipop_parameters = BIPOPParameters(
             self.lambda_, self.budget, self.mu / self.lambda_
         )
+        self.chiN = self.d ** 0.5 * (1 - 1 / (4 * self.d) + 1 / (21 * self.d ** 2))
+        self.ds = 2 - (2 / self.d)
+        self.beta = np.log(2) / max((np.sqrt(self.d) * np.log(self.d)), 1)
+        self.succes_ratio = .25
 
     def init_selection_parameters(self) -> None:
         """Initialization function for parameters that influence in selection."""
@@ -391,39 +398,39 @@ def init_adaptation_parameters(self) -> None:
 
         Examples are recombination weights and learning rates for the covariance
         matrix adapation.
-        TODO: clean the initlialization of these weights, this can be wrong in some cases
-        when mu != .5
         """
         if self.weights_option == "equal":
-            ws = np.ones(self.lambda_) / self.lambda_
-            self.weights = np.append(ws[:self.mu], ws[self.mu::-1] * -1)
-
-            if self.lambda_ % 2 != 0:
-                self.weights = np.append([1 / self.lambda_], self.weights[:-1])
+            ws = 1 / self.mu
+            self.weights = np.append(
+                np.ones(self.mu) * ws, np.ones(self.lambda_ - self.mu) * ws * -1
+            )
         elif self.weights_option == "1/2^lambda":
-            ws = 1 / 2 ** np.arange(1, self.lambda_ + 1) + (
-                (1 / (2 ** self.lambda_)) / self.lambda_
+            base = np.float64(2)
+            positive = self.mu / (base ** np.arange(1, self.mu + 1)) + (
+                (1 / (base ** self.mu)) / self.mu
             )
-            self.weights = np.append(ws[:self.mu], ws[self.mu::-1] * -1)
-            if self.lambda_ % 2 != 0:
-                self.weights = np.append([1 / self.lambda_ ** 2], self.weights[:-1])
+            n = self.lambda_ - self.mu
+            negative = (1 / (base ** np.arange(1, n + 1)) + (
+                (1 / (base ** n)) / n
+            ))[::-1] * -1
+            self.weights = np.append(positive, negative)
         else:
             self.weights = np.log((self.lambda_ + 1) / 2) - np.log(
                 np.arange(1, self.lambda_ + 1)
             )
 
         self.pweights = self.weights[: self.mu]
         self.nweights = self.weights[self.mu:]
-
         self.mueff = self.pweights.sum() ** 2 / (self.pweights ** 2).sum()
         mueff_neg = self.nweights.sum() ** 2 / (self.nweights ** 2).sum()
+
+        self.pweights = self.pweights / self.pweights.sum()
         self.c1 = self.c1 or 2 / ((self.d + 1.3) ** 2 + self.mueff)
         self.cmu = self.cmu or min(1 - self.c1, (2 * (
             (self.mueff - 2 + (1 / self.mueff))
             / ((self.d + 2) ** 2 + (2 * self.mueff / 2))
         )))
 
-        self.pweights = self.pweights / self.pweights.sum()
         amu_neg = 1 + (self.c1 / self.mu)
         amueff_neg = 1 + ((2 * mueff_neg) / (self.mueff + 2))
         aposdef_neg = (1 - self.c1 - self.cmu) / (self.d * self.cmu)
@@ -449,28 +456,26 @@ def init_adaptation_parameters(self) -> None:
         self.damps = 1.0 + (
             2.0 * max(0.0, np.sqrt((self.mueff - 1) / (self.d + 1)) - 1) + self.cs
         )
-        self.chiN = self.d ** 0.5 * (1 - 1 / (4 * self.d) + 1 / (21 * self.d ** 2))
-        self.ds = 2 - (2 / self.d)
-
-        self.beta = np.log(2) / max((np.sqrt(self.d) * np.log(self.d)), 1)
-        self.succes_ratio = .25
 
     def init_dynamic_parameters(self) -> None:
         """Initialization function of parameters that represent the dynamic state of the CMA-ES.
 
         Examples of such parameters are the Covariance matrix C and its 
         eigenvectors and the learning rate sigma.
         """
-        self.sigma = self.init_sigma
-        self.m = np.random.rand(self.d, 1)
-        self.m_old = np.empty((self.d, 1))
-        self.dm = np.zeros(self.d)
-        self.pc = np.zeros((self.d, 1))
-        self.ps = np.zeros((self.d, 1))
-        self.B = np.eye(self.d)
-        self.C = np.eye(self.d)
-        self.D = np.ones((self.d, 1))
-        self.inv_root_C = np.eye(self.d)
+        self.sigma = np.float64(self.sigma0)
+        if hasattr(self, "m") or self.x0 is None: 
+            self.m = np.float64(np.random.rand(self.d, 1))
+        else:
+            self.m = np.float64(self.x0.copy())
+        self.m_old = np.empty((self.d, 1), dtype=np.float64)
+        self.dm = np.zeros(self.d, dtype=np.float64)
+        self.pc = np.zeros((self.d, 1), dtype=np.float64)
+        self.ps = np.zeros((self.d, 1), dtype=np.float64)
+        self.B = np.eye(self.d, dtype=np.float64)
+        self.C = np.eye(self.d, dtype=np.float64)
+        self.D = np.ones((self.d, 1), dtype=np.float64)
+        self.inv_root_C = np.eye(self.d, dtype=np.float64)
         self.s = 0
         self.rank_tpa = None
         self.hs = True
@@ -594,11 +599,16 @@ def perform_eigendecomposition(self) -> None:
         ):
             self.init_dynamic_parameters()
         else:
+            C = self.C.copy()
             self.C = np.triu(self.C) + np.triu(self.C, 1).T
             self.D, self.B = linalg.eigh(self.C)
-            self.D = np.sqrt(self.D.astype(complex).reshape(-1, 1)).real
-            self.inv_root_C = np.dot(self.B, self.D ** -1 * self.B.T)
+            if np.all(self.D > 0):
+                self.D = np.sqrt(self.D.reshape(-1, 1))
+                self.inv_root_C = np.dot(self.B, self.D ** -1 * self.B.T)
+            else:
+                self.init_dynamic_parameters()
 
+
     def adapt_evolution_paths(self) -> None:
         """Method to adapt the evolution paths ps and pc."""
         self.dm = (self.m - self.m_old) / self.sigma
@@ -618,7 +628,7 @@ def adapt_evolution_paths(self) -> None:
     def perform_local_restart(self) -> None:
         """Method performing local restart, if a restart strategy is specified."""
         if self.local_restart:
-            if self.local_restart == "IPOP":
+            if self.local_restart == "IPOP" and self.mu > 512: 
                 self.mu *= self.ipop_factor
                 self.lambda_ *= self.ipop_factor
 
@@ -751,7 +761,7 @@ def calculate_termination_criteria(self) -> None:
         if self.local_restart or self.compute_termination_criteria:
             _t = self.t % self.d
             diag_C = np.diag(self.C.T)
-            d_sigma = self.sigma / self.init_sigma
+            d_sigma = self.sigma / self.sigma0
             best_fopts = self.best_fitnesses[self.last_restart:]
             median_fitnesses = self.median_fitnesses[self.last_restart:]
 
@@ -771,7 +781,7 @@ def calculate_termination_criteria(self) -> None:
                     ),
                     "tolx": np.all(
                         (np.append(self.pc.T, diag_C) * d_sigma)
-                        < (self.tolx * self.init_sigma)
+                        < (self.tolx * self.sigma0)
                     ),
                     "tolupsigma": (d_sigma > self.tolup_sigma * np.sqrt(self.D.max())),
                     "conditioncov": np.linalg.cond(self.C) > self.condition_cov,

modcma/utils.py

@@ -253,5 +253,5 @@ def ert(evals, budget):
             n_succ = (evals < budget).sum()
             _ert = float(evals.sum()) / int(n_succ)
         return _ert, np.std(evals), n_succ
-    except Exception:
+    except ZeroDivisionError:
         return float("inf"), np.nan, 0

tests/expected.py

@@ -17,14 +17,14 @@
         45.175034363167605,
         -52.01531942754774,
         1017.3030673363197,
-        79.03702012530019,
-        11.753343447706664,
-        131.38805502911111,
+        71.43871188992895,
+        3.146564668875037 ,
+        -8.279403543275532,
         -102.10176326857986,
         -543.6999174336394,
         42.700744146686965,
         -999.9997057338787,
-        8.695689533892846,
+        13.04988097465894,
         117.29640024093563,
     ],
     "base_sampler_gaussian": [
@@ -626,30 +626,30 @@
         117.58369741868964,
     ],
     "weights_option_1/2^lambda": [
-        79.5630351924981,
-        4629.283234379299,
-        -440.14992398127674,
-        -458.5322418694519,
-        17.401736662231777,
-        2171.3837540807126,
+        79.57458524535643,
+        6277.488135347779,
+        -444.31128872130535,
+        -458.9430906155403,
+        17.70768251894154,
+        2266.7358410247793,
         93.54091559010935,
-        149.3215688070224,
+        149.48997026626708,
         123.88066647543144,
-        4775.304616031011,
-        2485.875547646586,
-        2552406.7655545343,
-        57.15312947598912,
-        -51.89261437835274,
+        3477.138762872567,
+        3086.6572842078576,
+        2758225.2627536957,
+        59.257202521725354,
+        -51.89091502074251,
         1017.7040292314994,
-        73.89054993912737,
-        -12.20385658566134,
-        9.78378595249066,
-        -102.15555135313454,
+        72.22771406842722,
+        3.1834296917116482,
+        12.043970710043297,
+        -102.10653687894265,
         -543.698445859984,
         42.700744146686965,
-        -999.9821332012928,
-        18.573719726473,
-        111.84525007583429,
+        -999.9818092494717,
+        13.787357550858435,
+        111.45022572520682,
     ],
     "weights_option_default": [
         79.60197400323416,

tests/test_parameters.py

@@ -54,7 +54,7 @@ def test_bipop_parameters(self):
         self.p.used_budget += 11
         self.p.bipop_parameters.adapt(self.p.used_budget)
         self.assertEqual(self.p.bipop_parameters.large, True)
-        bp = self.p.bipop_parameters
+        bp = self.p.bipop_parameters 
         self.assertEqual(bp.lambda_, self.p.lambda_ * 2)
         self.assertEqual(bp.mu, self.p.mu * 2)
         self.assertEqual(bp.sigma, 2)
@@ -63,7 +63,7 @@ def test_bipop_parameters(self):
         self.assertEqual(self.p.bipop_parameters.large, False)
         self.assertLessEqual(bp.lambda_, self.p.lambda_)
         self.assertLessEqual(bp.mu, self.p.mu)
-        self.assertLessEqual(bp.sigma, self.p.init_sigma)
+        self.assertLessEqual(bp.sigma, self.p.sigma0)
         self.p.used_budget += 11
         bp.adapt(self.p.used_budget)
         self.assertEqual(bp.used_budget, 33)