implemented stch scalarisation

piglot/bin/piglot_plot.py

@@ -352,12 +352,16 @@ def plot_pareto(args):
             nondominated.append((point, variances[i, :] if has_variance else None))
         else:
             dominated.append((point, variances[i, :] if has_variance else None))
-    # Build the Pareto hull
-    hull = ConvexHull(pareto)
-    for simplex in hull.simplices:
-        # Hacky: filter the line between the two endpoints (assuming the list is sorted by x)
-        if abs(simplex[0] - simplex[1]) < pareto.shape[0] - 1:
-            ax.plot(pareto[simplex, 0], pareto[simplex, 1], 'r', ls='--')
+
+    # # Build the Pareto hull
+    # hull = ConvexHull(pareto)
+    # for simplex in hull.simplices:
+    #     # Hacky: filter the line between the two endpoints (assuming the list is sorted by x)
+    #     if abs(simplex[0] - simplex[1]) < pareto.shape[0] - 1:
+    #         ax.plot(pareto[simplex, 0], pareto[simplex, 1], 'b', ls='--', zorder=1000)
+
+    # Sort the points by x
+    nondominated = sorted(nondominated, key=lambda x: x[0][0])
     # Plot the points
     if has_variance:
         ax.errorbar(
@@ -369,6 +373,11 @@ def plot_pareto(args):
             fmt='o',
             label='Pareto front',
         )
+        ax.plot(
+            [point[0][0] for point in nondominated],
+            [point[0][1] for point in nondominated],
+            '--r',
+        )
         if args.all:
             ax.errorbar(
                 [point[0][0] for point in dominated],
@@ -385,13 +394,21 @@ def plot_pareto(args):
             [point[0][1] for point in nondominated],
             c='r',
             label='Pareto front',
+            zorder=1000,
+        )
+        ax.plot(
+            [point[0][0] for point in nondominated],
+            [point[0][1] for point in nondominated],
+            '--r',
+            zorder=1000,
         )
         if args.all:
             ax.scatter(
                 [point[0][0] for point in dominated],
                 [point[0][1] for point in dominated],
                 c='k',
                 label='Dominated points',
+                zorder=500,
             )
     if args.log:
         ax.set_xscale('log')

piglot/objective.py

@@ -150,7 +150,11 @@ class ObjectiveResult:
     """Container for objective results."""
     params: np.ndarray
     values: List[np.ndarray]
+    scalarisation: str
     variances: Optional[List[np.ndarray]] = None
+    weights: Optional[np.ndarray] = None
+    bounds: Optional[np.ndarray] = None
+    types: Optional[List[bool]] = None
 
     def __mc_variance(
         self,
@@ -180,6 +184,25 @@ def __mc_variance(
         ]
         return np.var(mc_objectives, axis=0)
 
+    @staticmethod
+    def normalise_objective(objective: np.ndarray, bounds: np.ndarray) -> np.ndarray:
+        """Normalise the objective.
+
+        Parameters
+        ----------
+        objective : np.ndarray
+            Objective to normalise.
+        bounds : np.ndarray
+            Bounds for the objectives.
+
+        Returns
+        -------
+        np.ndarray
+            Normalised objective.
+        """
+        # Normalise the objective
+        return (objective - bounds[:, 0]) / (bounds[:, 1] - bounds[:, 0])
+
     def scalarise(self, composition: Composition = None) -> float:
         """Scalarise the result under noise-free single-objective optimisation.
 
@@ -194,7 +217,32 @@ def scalarise(self, composition: Composition = None) -> float:
             Scalarised result.
         """
         if composition is None:
+            # Sanitise scalarisation method
+            if self.scalarisation not in ['mean', 'stch']:
+                raise ValueError(
+                    f"Invalid scalarisation '{self.scalarisation}'. Use 'mean' or 'stch'."
+                )
+
+            if self.scalarisation == "stch":
+                # Set all the objectives to be positive
+                values = abs(np.array(self.values))
+                # Set the weights, bounds and types
+                weights = np.array(self.weights)
+                bounds = np.array(self.bounds)
+                types = np.array(self.types)
+                # Calculate the costs and ideal point
+                costs = np.where(types, -1, 1)
+                ideal_point = np.where(types, 1, 0)
+                # Smoothing parameter for STCH
+                u = 0.01
+                # Calculate the normalised objective values
+                norm_funcs = self.normalise_objective(values, bounds)
+                # Calculate the Tchebycheff function value
+                tch_values = (np.abs((norm_funcs - ideal_point) * costs) / u) * weights
+                return np.log(np.sum(np.exp(tch_values))) * u
+
             return np.mean(self.values)
+
         return composition.composition(self.values, self.params).item()
 
     def scalarise_mo(self, composition: Composition = None) -> List[float]:

piglot/objectives/design.py

@@ -25,6 +25,7 @@ def __init__(
             negate: bool = False,
             weight: float = 1.0,
             n_points: int = None,
+            bounds: List[float] = None,
             ) -> None:
         # Sanitise prediction field
         if isinstance(prediction, str):
@@ -37,6 +38,8 @@ def __init__(
         self.weight = weight
         self.n_points = n_points
         self.flatten_utility = EndpointFlattenUtility(n_points) if n_points is not None else None
+        self.bounds = bounds
+        self.negate = negate
 
     @staticmethod
     def read(name: str, config: Dict[str, Any], output_dir: str) -> DesignTarget:
@@ -69,6 +72,7 @@ def read(name: str, config: Dict[str, Any], output_dir: str) -> DesignTarget:
             negate=bool(config.get('negate', False)),
             weight=float(config.get('weight', 1.0)),
             n_points=int(config['n_points']) if 'n_points' in config else None,
+            bounds=config.get('bounds', None),
         )
 
 
@@ -129,19 +133,26 @@ def __init__(
         composite: bool = False,
         multi_objective: bool = False,
         transformers: Dict[str, ResponseTransformer] = None,
+        scalarisation: str = 'mean',
     ) -> None:
         super().__init__(
             parameters,
             stochastic,
             composition=(
-                self.__composition(not multi_objective, stochastic, targets) if composite else None
+                self.__composition(
+                    not multi_objective,
+                    stochastic,
+                    targets,
+                    scalarisation
+                ) if composite else None
             ),
             num_objectives=len(targets) if multi_objective else 1,
             output_dir=output_dir,
         )
         self.solver = solver
         self.targets = targets
         self.transformers = transformers if transformers is not None else {}
+        self.scalarisation = scalarisation
 
     def prepare(self) -> None:
         """Prepare the objective for optimisation."""
@@ -153,6 +164,7 @@ def __composition(
         scalarise: bool,
         stochastic: bool,
         targets: List[DesignTarget],
+        scalarisation: str,
     ) -> Composition:
         """Build the composition utility for the design objective.
 
@@ -176,12 +188,19 @@ def __composition(
                 raise ValueError(
                     "All targets must have a number of points specified for the composition."
                 )
+        # Sanitise scalarisation method
+        if scalarisation not in ['mean', 'stch']:
+            raise ValueError(f"Invalid scalarisation '{scalarisation}'. Use 'mean' or 'stch'.")
+
         return ResponseComposition(
             scalarise=scalarise,
             stochastic=stochastic,
             weights=[t.weight for t in targets],
             reductions=[t.quantity for t in targets],
             flatten_list=[t.flatten_utility for t in targets],
+            scalarisation=scalarisation,
+            bounds=[t.bounds for t in targets] if scalarisation == 'stch' else None,
+            types=[t.negate for t in targets] if scalarisation == 'stch' else None,
         )
 
     @staticmethod
@@ -223,6 +242,10 @@ def _objective(self, values: np.ndarray, concurrent: bool = False) -> ObjectiveR
         ObjectiveResult
             Objective result.
         """
+        # Sanitise scalarisation method
+        if self.scalarisation not in ['mean', 'stch']:
+            raise ValueError(f"Invalid scalarisation '{self.scalarisation}'. Use 'mean' or 'stch'.")
+
         raw_responses = self.solver.solve(values, concurrent)
         # Transform responses
         for name, transformer in self.transformers.items():
@@ -252,7 +275,16 @@ def _objective(self, values: np.ndarray, concurrent: bool = False) -> ObjectiveR
             # Build statistical model for the target
             results.append(target.weight * np.mean(targets))
             variances.append(target.weight * np.var(targets) / len(targets))
-        return ObjectiveResult(values, results, variances if self.stochastic else None)
+
+        return ObjectiveResult(
+            values,
+            results,
+            self.scalarisation,
+            variances if self.stochastic else None,
+            [t.weight for t in self.targets] if self.scalarisation == 'stch' else None,
+            [t.bounds for t in self.targets] if self.scalarisation == 'stch' else None,
+            [t.negate for t in self.targets] if self.scalarisation == 'stch' else None,
+        )
 
     def plot_case(self, case_hash: str, options: Dict[str, Any] = None) -> List[Figure]:
         """Plot a given function call given the parameter hash
@@ -390,4 +422,5 @@ def read(
             composite=bool(config.get('composite', False)),
             multi_objective=bool(config.get('multi_objective', False)),
             transformers=transformers,
+            scalarisation=str(config.get('scalarisation', 'mean')),
         )

piglot/utils/composition/responses.py

@@ -295,6 +295,9 @@ def __init__(
         weights: List[float],
         reductions: List[Reduction],
         flatten_list: List[FlattenUtility],
+        scalarisation: str,
+        bounds: np.ndarray = None,
+        types: List[bool] = None,
     ) -> None:
         if len(flatten_list) != len(reductions):
             raise ValueError("Mismatched number of reductions and responses.")
@@ -305,6 +308,9 @@ def __init__(
         self.flatten_list = flatten_list
         self.lenghts = [flatten.length() for flatten in self.flatten_list]
         self.concat = ConcatUtility(self.lenghts)
+        self.bounds = bounds
+        self.scalarisation = scalarisation
+        self.types = types
 
     def transform(self, params: np.ndarray, responses: List[List[OutputResult]]) -> ObjectiveResult:
         """Transform a set of responses into a fixed-size ObjectiveResult for the optimiser.
@@ -364,6 +370,27 @@ def _expand_params(time: torch.Tensor, params: torch.Tensor) -> torch.Tensor:
         # Expand the parameters along the first dimensions
         return params.expand(*([a for a in time.shape[:-1]] + [params.shape[-1]]))
 
+    @staticmethod
+    def normalise_objective(objective: torch.Tensor, bounds: torch.Tensor) -> torch.Tensor:
+        """Normalise the objective.
+
+        Parameters
+        ----------
+        objective : torch.Tensor
+            Objective to normalise.
+        bounds : torch.Tensor
+            Bounds for the objectives.
+        types : torch.Tensor
+            Types for the objectives.
+
+        Returns
+        -------
+        torch.Tensor
+            Normalised objective.
+        """
+        # Normalise the objective
+        return (objective - bounds[:, 0]) / (bounds[:, 1] - bounds[:, 0])
+
     def composition_torch(self, inner: torch.Tensor, params: torch.Tensor) -> torch.Tensor:
         """Compute the outer function of the composition with gradients.
 
@@ -379,6 +406,10 @@ def composition_torch(self, inner: torch.Tensor, params: torch.Tensor) -> torch.
         torch.Tensor
             Composition result.
         """
+        # Sanitise the scalarisation method
+        if self.scalarise and self.scalarisation not in ['mean', 'stch']:
+            raise ValueError(f"Invalid scalarisation '{self.scalarisation}'. Use 'mean' or 'stch'.")
+
         # Split the inner responses
         responses = self.concat.split_torch(inner)
         # Unflatten each response
@@ -391,7 +422,30 @@ def composition_torch(self, inner: torch.Tensor, params: torch.Tensor) -> torch.
             reduction.reduce_torch(time, data, self._expand_params(time, params))
             for (time, data), reduction in zip(unflattened, self.reductions)
         ], dim=-1)
-        # Apply the weights
-        objective = objective * torch.tensor(self.weights).to(inner.device)
-        # If needed, scalarise the objectives
-        return torch.mean(objective, dim=-1) if self.scalarise else objective
+
+        # Mean scalarisation if requested
+        if self.scalarise:
+            # Smooth TCHebycheff scalarisation (STCH) if requested
+            if self.scalarisation == 'stch':
+                # Set all the objectives to be positive
+                objective = objective.abs()
+                # Set the weights, bounds and types
+                weights = torch.tensor(self.weights).to(inner.device)
+                bounds = torch.tensor(self.bounds).to(inner.device)
+                types = torch.tensor(self.types).to(inner.device)
+                # Calculate the costs and ideal point
+                costs = torch.where(types, torch.tensor(-1.0), torch.tensor(1.0))
+                ideal_point = torch.where(types, torch.tensor(1.0), torch.tensor(0.0))
+                # Smoothing parameter for STCH
+                u = 0.01
+                # Calculate the normalised objective values
+                norm_objective = self.normalise_objective(objective, bounds)
+                # Calculate the Tchebycheff function value
+                tch_values = (torch.abs((norm_objective - ideal_point) * costs) / u) * weights
+                return torch.log(torch.sum(torch.exp(tch_values), dim=-1)) * u
+            # Mean scalarisation otherwise
+            # Apply the weights
+            objective = objective * torch.tensor(self.weights).to(inner.device)
+            return torch.mean(objective, dim=-1)
+
+        return objective * torch.tensor(self.weights).to(inner.device)