diff --git a/pybop/_optimisation.py b/pybop/_optimisation.py index 67a40a151..63821735d 100644 --- a/pybop/_optimisation.py +++ b/pybop/_optimisation.py @@ -22,7 +22,7 @@ class Optimisation: If True, the optimization progress is printed (default: False). physical_viability : bool, optional If True, the feasibility of the optimised parameters is checked (default: True). - allow_infeasible_solutions : bool, optional + allow_infeasible_solutions : bool or list[bool], optional If True, infeasible parameter values will be allowed in the optimisation (default: True). Attributes @@ -65,9 +65,14 @@ def __init__( # Set whether to allow infeasible locations if self.cost.problem is not None and hasattr(self.cost.problem, "_model"): - self.cost.problem._model.allow_infeasible_solutions = ( - self.allow_infeasible_solutions - ) + try: + for model, allow_infeasible in zip( + self.cost.problem._model, self.allow_infeasible_solutions + ): + model.allow_infeasible_solutions = allow_infeasible + except TypeError: + for model in self.cost.problem._model: + model.allow_infeasible_solutions = self.allow_infeasible_solutions else: # Turn off this feature as there is no model self.physical_viability = False @@ -517,13 +522,16 @@ def check_optimal_parameters(self, x): Check if the optimised parameters are physically viable. """ - if self.cost.problem._model.check_params( - inputs=x, allow_infeasible_solutions=False + if all( + [ + model.check_params(inputs=x, allow_infeasible_solutions=False) + for model in self.cost.problem._model + ] ): return else: warnings.warn( - "Optimised parameters are not physically viable! \nConsider retrying the optimisation" + "One or more optimised parameters are not physically viable! \nConsider retrying the optimisation" + " with a non-gradient-based optimiser and the option allow_infeasible_solutions=False", UserWarning, stacklevel=2, diff --git a/pybop/costs/fitting_costs.py b/pybop/costs/fitting_costs.py index 0665454b0..4b9ac5d5c 100644 --- a/pybop/costs/fitting_costs.py +++ b/pybop/costs/fitting_costs.py @@ -17,8 +17,8 @@ class RootMeanSquaredError(BaseCost): """ - def __init__(self, problem): - super(RootMeanSquaredError, self).__init__(problem) + def __init__(self, problem, weights=None): + super(RootMeanSquaredError, self).__init__(problem, weights=weights) # Default fail gradient self._de = 1.0 @@ -41,62 +41,61 @@ def _evaluate(self, x, grad=None): The root mean square error. """ - prediction = self.problem.evaluate(x) - - for key in self.signal: - if len(prediction.get(key, [])) != len(self._target.get(key, [])): - return np.float64(np.inf) # prediction doesn't match target - - e = np.array( - [ - np.sqrt(np.mean((prediction[signal] - self._target[signal]) ** 2)) - for signal in self.signal - ] - ) + predictions = self.problem.evaluate(x) + errors = np.zeros(len(predictions)) + # For each dataset... + for i, (prediction, target) in enumerate(zip(predictions, self._target)): + for key in self.signal: + if len(prediction.get(key, [])) != len(target.get(key, [])): + return np.float64(np.inf) # prediction doesn't match target + + # Calculate the root mean square error + e = np.array( + [ + np.sqrt(np.mean((prediction[signal] - target[signal]) ** 2)) + for signal in self.signal + ] + ) + if self.n_outputs == 1: + errors[i] = e.item() + else: + errors[i] = np.sum(e) - if self.n_outputs == 1: - return e.item() - else: - return np.sum(e) + # Weight errors and return + if self.weights is None: + return errors.sum() + return np.inner(self.weights, errors) def _evaluateS1(self, x): - """ - Compute the cost and its gradient with respect to the parameters. - - Parameters - ---------- - x : array-like - The parameters for which to compute the cost and gradient. + ys, dys = self.problem.evaluateS1(x) + + errors, gradients = [], [] + for i, (prediction, dy, target) in enumerate(zip(ys, dys, self._target)): + for key in self.signal: + if len(prediction.get(key, [])) != len(target.get(key, [])): + e = np.float64(np.inf) + de = self._de * np.ones(self.n_parameters) + return e, de + + r = np.array( + [prediction[signal] - target[signal] for signal in self.signal] + ) + errors.append(np.sqrt(np.mean(r**2, axis=1))) + gradients.append( + np.mean((r * dy.T), axis=2) + / (np.sqrt(np.mean((r * dy.T) ** 2, axis=2)) + np.finfo(float).eps) + ) - Returns - ------- - tuple - A tuple containing the cost and the gradient. The cost is a float, - and the gradient is an array-like of the same length as `x`. + errors = np.array(errors) + gradients = np.array(gradients) + if self.weights is None: + return errors.sum(), gradients.sum(axis=1) - Raises - ------ - ValueError - If an error occurs during the calculation of the cost or gradient. - """ - y, dy = self.problem.evaluateS1(x) - - for key in self.signal: - if len(y.get(key, [])) != len(self._target.get(key, [])): - e = np.float64(np.inf) - de = self._de * np.ones(self.n_parameters) - return e, de - - r = np.array([y[signal] - self._target[signal] for signal in self.signal]) - e = np.sqrt(np.mean(r**2, axis=1)) - de = np.mean((r * dy.T), axis=2) / ( - np.sqrt(np.mean((r * dy.T) ** 2, axis=2)) + np.finfo(float).eps + # Weight and return + weighted_gradients = np.inner(self.weights.reshape((1, 1, -1)), gradients).sum( + axis=1 ) - - if self.n_outputs == 1: - return e.item(), de.flatten() - else: - return np.sum(e), np.sum(de, axis=1) + return np.inner(self.weights, errors), weighted_gradients def set_fail_gradient(self, de): """ @@ -155,22 +154,29 @@ def _evaluate(self, x, grad=None): float The sum of squared errors. """ - prediction = self.problem.evaluate(x) - - for key in self.signal: - if len(prediction.get(key, [])) != len(self._target.get(key, [])): - return np.float64(np.inf) # prediction doesn't match target + predictions = self.problem.evaluate(x) + errors = np.zeros(len(predictions)) + # For each dataset... + for i, (prediction, target) in enumerate(zip(predictions, self._target)): + for key in self.signal: + if len(prediction.get(key, [])) != len(target.get(key, [])): + return np.float64(np.inf) # prediction doesn't match target + + e = np.array( + [ + np.sum(((prediction[signal] - target[signal]) ** 2)) + for signal in self.signal + ] + ) + if self.n_outputs == 1: + errors[i] = e.item() + else: + errors[i] = np.sum(e) - e = np.array( - [ - np.sum(((prediction[signal] - self._target[signal]) ** 2)) - for signal in self.signal - ] - ) - if self.n_outputs == 1: - return e.item() - else: - return np.sum(e) + # Weight errors and return + if self.weights is None: + return errors.sum() + return np.inner(self.weights, errors) def _evaluateS1(self, x): """ @@ -192,18 +198,29 @@ def _evaluateS1(self, x): ValueError If an error occurs during the calculation of the cost or gradient. """ - y, dy = self.problem.evaluateS1(x) - for key in self.signal: - if len(y.get(key, [])) != len(self._target.get(key, [])): - e = np.float64(np.inf) - de = self._de * np.ones(self.n_parameters) - return e, de - - r = np.array([y[signal] - self._target[signal] for signal in self.signal]) - e = np.sum(np.sum(r**2, axis=0), axis=0) - de = 2 * np.sum(np.sum((r * dy.T), axis=2), axis=1) - - return e, de + ys, dys = self.problem.evaluateS1(x) + + errors, gradients = [], [] + for i, (prediction, dy, target) in enumerate(zip(ys, dys, self._target)): + for key in self.signal: + if len(prediction.get(key, [])) != len(target.get(key, [])): + e = np.float64(np.inf) + de = self._de * np.ones(self.n_parameters) + return e, de + + r = np.array( + [prediction[signal] - target[signal] for signal in self.signal] + ) + errors.append(np.sum(np.sum(r**2, axis=0), axis=0)) + gradients.append(2 * np.sum(np.sum((r * dy.T), axis=2), axis=1)) + + errors = np.array(errors) + gradients = np.array(gradients) + if self.weights is None: + return errors, gradients + return np.inner(self.weights, errors), np.inner(self.weights, gradients) + + def set_fail_gradient(self, de): """ diff --git a/pybop/plotting/plot_dataset.py b/pybop/plotting/plot_dataset.py index 483d7035a..a4a9952a8 100644 --- a/pybop/plotting/plot_dataset.py +++ b/pybop/plotting/plot_dataset.py @@ -33,30 +33,34 @@ def plot_dataset( # Get data dictionary dataset.check(signal) - # Compile ydata and labels or legend - y = [dataset[s] for s in signal] - if len(signal) == 1: - yaxis_title = signal[0] - if trace_names is None: - trace_names = ["Data"] - else: - yaxis_title = "Output" - if trace_names is None: - trace_names = pybop.StandardPlot.remove_brackets(signal) - - # Create the figure - fig = pybop.plot_trajectories( - x=dataset["Time [s]"], - y=y, - trace_names=trace_names, - show=False, - xaxis_title="Time / s", - yaxis_title=yaxis_title, - ) - fig.update_layout(**layout_kwargs) - if "ipykernel" in sys.modules and show: - fig.show("svg") - elif show: - fig.show() + if not hasattr(dataset, '__len__'): + dataset = [dataset] + + for data in dataset: + # Compile ydata and labels or legend + y = [data[s] for s in signal] + if len(signal) == 1: + yaxis_title = signal[0] + if trace_names is None: + trace_names = ["Data"] + else: + yaxis_title = "Output" + if trace_names is None: + trace_names = pybop.StandardPlot.remove_brackets(signal) + + # Create the figure + fig = pybop.plot_trajectories( + x=data["Time [s]"], + y=y, + trace_names=trace_names, + show=False, + xaxis_title="Time / s", + yaxis_title=yaxis_title, + ) + fig.update_layout(**layout_kwargs) + if "ipykernel" in sys.modules and show: + fig.show("svg") + elif show: + fig.show() return fig diff --git a/pybop/plotting/plot_problem.py b/pybop/plotting/plot_problem.py index 8540eb4c5..b65380192 100644 --- a/pybop/plotting/plot_problem.py +++ b/pybop/plotting/plot_problem.py @@ -34,62 +34,65 @@ def quick_plot(problem, parameter_values=None, show=True, **layout_kwargs): parameter_values = problem.x0 # Extract the time data and evaluate the model's output and target values - xaxis_data = problem.time_data() - model_output = problem.evaluate(parameter_values) - target_output = problem.get_target() + xaxis_datasets = problem.time_data() + model_outputs = problem.evaluate(parameter_values) + target_outputs = problem.get_target() # Create a plot for each output figure_list = [] for i in problem.signal: - default_layout_options = dict( - title="Scatter Plot", - xaxis_title="Time / s", - yaxis_title=pybop.StandardPlot.remove_brackets(i), - ) - - # Create a plotting dictionary - if isinstance(problem, pybop.DesignProblem): - trace_name = "Optimised" - opt_time_data = model_output["Time [s]"] - else: - trace_name = "Model" - opt_time_data = xaxis_data - - plot_dict = pybop.StandardPlot( - x=opt_time_data, - y=model_output[i], - layout_options=default_layout_options, - trace_names=trace_name, - ) - - target_trace = plot_dict.create_trace( - x=xaxis_data, - y=target_output[i], - name="Reference", - mode="markers", - showlegend=True, - ) - plot_dict.traces.append(target_trace) - - if isinstance(problem, pybop.FittingProblem): - # Compute the standard deviation as proxy for uncertainty - plot_dict.sigma = np.std(model_output[i] - target_output[i]) - - # Convert x and upper and lower limits into lists to create a filled trace - x = xaxis_data.tolist() - y_upper = (model_output[i] + plot_dict.sigma).tolist() - y_lower = (model_output[i] - plot_dict.sigma).tolist() - - fill_trace = plot_dict.create_trace( - x=x + x[::-1], - y=y_upper + y_lower[::-1], - fill="toself", - fillcolor="rgba(255,229,204,0.8)", - line=dict(color="rgba(255,255,255,0)"), - hoverinfo="skip", - showlegend=False, + for target_output, model_output, xaxis_data in zip( + target_outputs, model_outputs, xaxis_datasets + ): + default_layout_options = dict( + title="Scatter Plot", + xaxis_title="Time / s", + yaxis_title=pybop.StandardPlot.remove_brackets(i), ) - plot_dict.traces.append(fill_trace) + + # Create a plotting dictionary + if isinstance(problem, pybop.DesignProblem): + trace_name = "Optimised" + opt_time_data = model_output["Time [s]"] + else: + trace_name = "Model" + opt_time_data = xaxis_data + + plot_dict = pybop.StandardPlot( + x=opt_time_data, + y=model_output[i], + layout_options=default_layout_options, + trace_names=trace_name, + ) + + target_trace = plot_dict.create_trace( + x=xaxis_data, + y=target_output[i], + name="Reference", + mode="markers", + showlegend=True, + ) + plot_dict.traces.append(target_trace) + + if isinstance(problem, pybop.FittingProblem): + # Compute the standard deviation as proxy for uncertainty + plot_dict.sigma = np.std(model_output[i] - target_output[i]) + + # Convert x and upper and lower limits into lists to create a filled trace + x = xaxis_data.tolist() + y_upper = (model_output[i] + plot_dict.sigma).tolist() + y_lower = (model_output[i] - plot_dict.sigma).tolist() + + fill_trace = plot_dict.create_trace( + x=x + x[::-1], + y=y_upper + y_lower[::-1], + fill="toself", + fillcolor="rgba(255,229,204,0.8)", + line=dict(color="rgba(255,255,255,0)"), + hoverinfo="skip", + showlegend=False, + ) + plot_dict.traces.append(fill_trace) # Reverse the order of the traces to put the model on top plot_dict.traces = plot_dict.traces[::-1]