diff --git a/piglot/bin/piglot_plot.py b/piglot/bin/piglot_plot.py index a9aa09d..3cdad6a 100644 --- a/piglot/bin/piglot_plot.py +++ b/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,6 +394,13 @@ 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( @@ -392,6 +408,7 @@ def plot_pareto(args): [point[0][1] for point in dominated], c='k', label='Dominated points', + zorder=500, ) if args.log: ax.set_xscale('log') diff --git a/piglot/objective.py b/piglot/objective.py index 29fee19..4aa44b7 100644 --- a/piglot/objective.py +++ b/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]: diff --git a/piglot/objectives/design.py b/piglot/objectives/design.py index 90cdc2c..c811a9b 100644 --- a/piglot/objectives/design.py +++ b/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,12 +133,18 @@ 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, @@ -142,6 +152,7 @@ def __init__( 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')), ) diff --git a/piglot/utils/composition/responses.py b/piglot/utils/composition/responses.py index ee32bf4..fa7ee8c 100644 --- a/piglot/utils/composition/responses.py +++ b/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)