Change Constraints Convention

src/paretobench/cf.py

@@ -64,7 +64,7 @@ def _call(self, x):
             )
         )
         g = np.vstack((f[0] + f[1] - self.a * np.abs(np.sin(self.b * np.pi * (f[0] - f[1] + 1))) - 1,))
-        return Population(f=f.T, g=g.T)
+        return Population(f=f.T, g=-g.T)
 
     @property
     def var_lower_bounds(self):
@@ -110,7 +110,7 @@ def _call(self, x):
         )
         t = f[1] + np.sqrt(f[0]) - self.a * np.sin(self.b * np.pi * (np.sqrt(f[0]) - f[1] + 1)) - 1
         g = np.vstack((t / (1 + np.exp(4 * np.abs(t))),))
-        return Population(f=f.T, g=g.T)
+        return Population(f=f.T, g=-g.T)
 
     @property
     def var_lower_bounds(self):
@@ -165,7 +165,7 @@ def _call(self, x):
             )
         )
         g = np.vstack((f[1] + f[0] ** 2 - self.a * np.sin(self.b * np.pi * (f[0] ** 2 - f[1] + 1)) - 1,))
-        return Population(f=f.T, g=g.T)
+        return Population(f=f.T, g=-g.T)
 
     @property
     def var_lower_bounds(self):
@@ -218,7 +218,7 @@ def _call(self, x):
 
         t = x[1] - np.sin(6 * np.pi * x[0] + 2 * np.pi / self.n) - 0.5 * x[0] + 0.25
         g = np.vstack((t / (1 + np.exp(4 * np.abs(t))),))
-        return Population(f=f.T, g=g.T)
+        return Population(f=f.T, g=-g.T)
 
     @property
     def var_lower_bounds(self):
@@ -270,7 +270,7 @@ def _call(self, x):
 
         g = np.vstack((x[1] - 0.8 * x[0] * np.sin(6 * np.pi * x[0] + 2 * np.pi / self.n) - 0.5 * x[0] + 0.25,))
 
-        return Population(f=f.T, g=g.T)
+        return Population(f=f.T, g=-g.T)
 
     @property
     def var_lower_bounds(self):
@@ -328,7 +328,7 @@ def _call(self, x):
         )
         g = np.vstack((g1, g2))
 
-        return Population(f=f.T, g=g.T)
+        return Population(f=f.T, g=-g.T)
 
     @property
     def var_lower_bounds(self):
@@ -384,7 +384,7 @@ def _call(self, x):
         )
         g = np.vstack((g1, g2))
 
-        return Population(f=f.T, g=g.T)
+        return Population(f=f.T, g=-g.T)
 
     @property
     def var_lower_bounds(self):
@@ -443,7 +443,7 @@ def _call(self, x):
             )
         )
 
-        return Population(f=f.T, g=g.T)
+        return Population(f=f.T, g=-g.T)
 
     @property
     def var_lower_bounds(self):
@@ -501,7 +501,7 @@ def _call(self, x):
             )
         )
 
-        return Population(f=f.T, g=g.T)
+        return Population(f=f.T, g=-g.T)
 
     @property
     def var_lower_bounds(self):
@@ -556,7 +556,7 @@ def _call(self, x):
             )
         )
 
-        return Population(f=f.T, g=g.T)
+        return Population(f=f.T, g=-g.T)
 
     @property
     def var_lower_bounds(self):

src/paretobench/containers.py

@@ -76,7 +76,7 @@ def set_default_vals(cls, values):
         if values.get("obj_directions") is None:
             values["obj_directions"] = "-" * values["f"].shape[1]
         if values.get("constraint_directions") is None:
-            values["constraint_directions"] = ">" * values["g"].shape[1]
+            values["constraint_directions"] = "<" * values["g"].shape[1]
         if values.get("constraint_targets") is None:
             values["constraint_targets"] = np.zeros(values["g"].shape[1], dtype=float)
 
@@ -199,8 +199,8 @@ def g_canonical(self):
         """
         Return constraints transformed such that g[...] >= 0 are the feasible solutions.
         """
-        gc = binary_str_to_numpy(self.constraint_directions, ">", "<")[None, :] * self.g
-        gc += binary_str_to_numpy(self.constraint_directions, "<", ">")[None, :] * self.constraint_targets[None, :]
+        gc = binary_str_to_numpy(self.constraint_directions, "<", ">")[None, :] * self.g
+        gc += binary_str_to_numpy(self.constraint_directions, ">", "<")[None, :] * self.constraint_targets[None, :]
         return gc
 
     def __add__(self, other: "Population") -> "Population":
@@ -289,7 +289,7 @@ def get_nondominated_indices(self):
     def get_feasible_indices(self):
         if self.g.shape[1] == 0:
             return np.ones((len(self)), dtype=bool)
-        return np.all(self.g_canonical >= 0.0, axis=1)
+        return np.all(self.g_canonical <= 0.0, axis=1)
 
     def get_nondominated_set(self):
         return self[self.get_nondominated_indices()]
@@ -619,7 +619,7 @@ def _to_h5py_group(self, g: h5py.Group):
             g["g"].attrs["targets"] = self.reports[0].constraint_targets
 
     @classmethod
-    def _from_h5py_group(cls, grp: h5py.Group):
+    def _from_h5py_group(cls, grp: h5py.Group, file_version: str):
         """
         Construct a new History object from data in an HDF5 group.
 
@@ -638,12 +638,17 @@ def _from_h5py_group(cls, grp: h5py.Group):
         f = grp["f"][()]
         g = grp["g"][()]
 
-        # Get the names
+        # Get the objective / constraint settings
         obj_directions = grp["f"].attrs.get("directions", None)
         constraint_directions = grp["g"].attrs.get("directions", None)
         constraint_targets = grp["g"].attrs.get("targets", None)
 
-        # Get the objective / constraint settings
+        # Before file version 1.1.0 which introduced explicit constraint directions,
+        # the default constraint type was g(x) >= 0.0
+        if file_version == "1.0.0":
+            constraint_directions = ">" * g.shape[1]
+
+        # Get the names
         names_x = grp["x"].attrs.get("names", None)
         names_f = grp["f"].attrs.get("names", None)
         names_g = grp["g"].attrs.get("names", None)
@@ -722,6 +727,10 @@ class Experiment(BaseModel):
     """
     Represents on "experiment" performed on a multibojective genetic algorithm. It may contain several evaluations of the
     algorithm on different problems or repeated iterations on the same problem.
+
+    Note: when loading files from disk, the `file_version` field will indicate the version of the file
+    that was loaded. When the `Experiment` is saved, the resulting file will contain the file version
+    used to save the data.
     """
 
     runs: List[History]
@@ -850,7 +859,7 @@ def save(self, fname):
             f.attrs["software_version"] = self.software_version
             f.attrs["comment"] = self.comment
             f.attrs["creation_time"] = self.creation_time.isoformat()
-            f.attrs["file_version"] = self.file_version
+            f.attrs["file_version"] = "1.1.0"
             f.attrs["file_format"] = "ParetoBench Multi-Objective Optimization Data"
 
             # Calculate the necessary zero padding based on the number of runs
@@ -886,11 +895,12 @@ def load(cls, fname):
         # Load the data
         with h5py.File(fname, mode="r") as f:
             # Load each of the runs keeping track of the order of the indices
+            file_version = f.attrs["file_version"]
             idx_runs = []
             for idx_str, run_grp in f.items():
                 m = re.match(r"run_(\d+)", idx_str)
                 if m:
-                    idx_runs.append((int(m.group(1)), History._from_h5py_group(run_grp)))
+                    idx_runs.append((int(m.group(1)), History._from_h5py_group(run_grp, file_version=file_version)))
             runs = [x[1] for x in sorted(idx_runs, key=lambda x: x[0])]
 
             # Convert the creation_time back to a timezone-aware datetime object
@@ -905,4 +915,5 @@ def load(cls, fname):
                 software_version=f.attrs["software_version"],
                 comment=f.attrs["comment"],
                 creation_time=creation_time,
+                file_version=file_version,
             )

src/paretobench/ctp.py

@@ -62,7 +62,7 @@ def _call(self, x):
         g = []
         for i in range(self.j):
             g.append(f[1] - self._a[i] * np.exp(-self._b[i] * f[0]))
-        return Population(f=f.T, g=np.vstack(g).T)
+        return Population(f=f.T, g=-np.vstack(g).T)
 
 
 class CTP2_7(CTPx):
@@ -118,7 +118,7 @@ def _call(self, x):
         )
         g = np.vstack((c,))
 
-        return Population(f=f.T, g=g.T)
+        return Population(f=f.T, g=-g.T)
 
 
 class CTP2(CTP2_7):

src/paretobench/dtlz.py

@@ -182,7 +182,7 @@ def _call(self, x):
         fsum = f[:-1, None, :] + f[None, :-1, :]
         fsum[np.diag_indices(2), :] = np.inf  # For i != j
         g = np.vstack(g + [2 * f[-1] + np.min(np.min(fsum, axis=0), axis=0) - 1])
-        return Population(f=f.T, g=g.T)
+        return Population(f=f.T, g=-g.T)
 
     def get_pareto_front(self, n):
         # Break points into the "pole" feature and the lower PF. Based on dimension, number in
@@ -244,7 +244,7 @@ def _call(self, x):
             ]
         )
         g = np.vstack([f[-1] ** 2 + f[j] ** 2 - 1 for j in range(self.m - 1)])
-        return Population(f=f.T, g=g.T)
+        return Population(f=f.T, g=-g.T)
 
     def get_pareto_front(self, n):
         th = np.linspace(0, 1, n)

src/paretobench/misc.py

@@ -196,7 +196,7 @@ def _call(self, x):
 
         f = np.array([x[0], (1 + x[1]) / x[0]])
         g = np.array([x[1] + 9 * x[0] - 6, -x[1] + 9 * x[0] - 1])
-        return Population(f=f.T, g=g.T)
+        return Population(f=f.T, g=-g.T)
 
     @property
     def var_lower_bounds(self):
@@ -233,7 +233,7 @@ def _call(self, x):
 
         f = np.array([(x[0] - 2) ** 2 + (x[1] - 1) ** 2 + 2, 9 * x[0] - (x[1] - 1) ** 2])
         g = np.array([225 - (x[0] ** 2 + x[1] ** 2), -10 - (x[0] - 3 * x[1])])
-        return Population(f=f.T, g=g.T)
+        return Population(f=f.T, g=-g.T)
 
     @property
     def var_lower_bounds(self):
@@ -275,7 +275,7 @@ def _call(self, x):
                 0.5 - ((x[0] - 0.5) ** 2 + (x[1] - 0.5) ** 2),
             ]
         )
-        return Population(f=f.T, g=g.T)
+        return Population(f=f.T, g=-g.T)
 
     @property
     def var_bounds(self):
@@ -334,7 +334,7 @@ def _call(self, x):
                 550.0 - (0.164 / (x[0] * x[1]) + 631.13 * x[2] - 54.48),
             ]
         )
-        return Population(f=f.T, g=g.T)
+        return Population(f=f.T, g=-g.T)
 
     @property
     def var_lower_bounds(self):

src/paretobench/utils.py

@@ -41,7 +41,23 @@ def uniform_grid(n, m):
     )
 
 
-def get_domination(objs, constraints=None):
+def get_domination(objs: np.ndarray, constraints: np.ndarray | None = None) -> np.ndarray:
+    """
+    Calculate domination between all pairs of individuals in the population.
+
+    Parameters
+    ----------
+    objs : np.ndarray
+        The objectives, first index is batch index
+    constraints : np.ndarray / None
+        The constraints (such that g<=0 means feasible) or None if no constraints
+
+    Returns
+    -------
+    np.ndarray
+        (n, n) boolean array with domination relationship. dom_{ij} is true if individual i
+        dominates individual j.
+    """
     # Compare all pairs of individuals based on domination
     dom = np.bitwise_and(
         (objs[:, None, :] <= objs[None, :, :]).all(axis=-1),
@@ -50,14 +66,14 @@ def get_domination(objs, constraints=None):
 
     if constraints is not None:
         # If one individual is feasible and the other isn't, set domination
-        feas = constraints >= 0.0
+        feas = constraints <= 0.0
         ind = np.bitwise_and(feas.all(axis=1)[:, None], ~feas.all(axis=1)[None, :])
         dom[ind] = True
         ind = np.bitwise_and(~feas.all(axis=1)[:, None], feas.all(axis=1)[None, :])
         dom[ind] = False
 
         # If both are infeasible, then the individual with the least constraint violation wins
-        constraint_violation = -np.sum(np.minimum(constraints, 0), axis=1)
+        constraint_violation = np.sum(np.maximum(constraints, 0), axis=1)
         comp = constraint_violation[:, None] < constraint_violation[None, :]
         ind = ~np.bitwise_or(feas.all(axis=1)[:, None], feas.all(axis=1)[None, :])
         dom[ind] = comp[ind]

tests/test_containers.py

@@ -94,20 +94,20 @@ def test_generate_population():
         (np.array([[1.0, 0.0], [0.0, 1.0]]).T, None, [True, True]),
         (np.array([[1.0, 0.0], [1.0, 1.0]]).T, None, [False, True]),
         (np.array([[0.0, 1.0], [1.0, 1.0]]).T, None, [True, False]),
-        (np.array([[0.0, 1.0], [1.0, 0.0]]).T, np.array([[1.0, 0.1]]).T, [True, True]),
+        (np.array([[0.0, 1.0], [1.0, 0.0]]).T, -np.array([[1.0, 0.1]]).T, [True, True]),
         (
             np.array([[0.0, 1.0], [1.0, 0.0]]).T,
-            np.array([[1.0, -0.1]]).T,
+            -np.array([[1.0, -0.1]]).T,
             [True, False],
         ),
         (
             np.array([[0.0, 1.0], [1.0, 0.0]]).T,
-            np.array([[-1.0, -0.1]]).T,
+            -np.array([[-1.0, -0.1]]).T,
             [False, True],
         ),
         (
             np.array([[0.0, 1.0], [1.0, 0.0]]).T,
-            np.array([[-1.0, -0.1], [100.0, -0.1], [100.0, 3e3]]).T,
+            -np.array([[-1.0, -0.1], [100.0, -0.1], [100.0, 3e3]]).T,
             [False, True],
         ),
     ],

tests/test_plotting.py

@@ -19,7 +19,7 @@ def test_per_point_settings_nondominated():
     # Create a simple population where all points are feasible and non-dominated
     pop = Population.from_random(n_objectives=2, n_decision_vars=2, n_constraints=1, pop_size=3)
     # Make all points feasible
-    pop.g[:] = 1.0
+    pop.g[:] = -1.0
     # Make points non-dominated by setting objectives such that no point dominates another
     pop.f = np.array([[1.0, 2.0], [2.0, 1.0], [1.5, 1.5]])
 
@@ -39,7 +39,7 @@ def test_per_point_settings_random():
     # Set up specific test cases
     pop.f = np.array([[1.0, 1.0], [2.0, 2.0], [1.5, 1.5], [0.5, 0.5], [1.5, 0.5]])
 
-    pop.g = np.array([[1.0], [1.0], [-1.0], [-1.0], [1.0]])
+    pop.g = -np.array([[1.0], [1.0], [-1.0], [-1.0], [1.0]])
 
     settings = get_per_point_settings_population(pop, domination_filt="all", feasibility_filt="all")