from __future__ import annotations
#
# Solution class
#
import json
import numbers
import pickle
from functools import cached_property
from itertools import chain
import casadi
import numpy as np
import pandas as pd
from scipy.io import savemat
import pybamm
from pybamm.codegen.compilation import aot_compile
class NumpyEncoder(json.JSONEncoder):
"""
Numpy serialiser helper class that converts numpy arrays to a list.
Numpy arrays cannot be directly converted to JSON, so the arrays are
converted to python list objects before encoding.
"""
def default(self, obj):
if isinstance(obj, np.ndarray):
return obj.tolist()
# won't be called since we only need to convert numpy arrays
return json.JSONEncoder.default(self, obj) # pragma: no cover
class SolutionBase:
"""Base class for all PyBaMM solution types (time-series and EIS).
Subclasses populate ``self._data`` (a :class:`~pybamm.FuzzyDict`) with
named data arrays so that ``solution["Variable name"]`` works uniformly
across solution types.
"""
def __init__(self):
self.set_up_time = None
self.solve_time = None
self._data = pybamm.FuzzyDict()
def __getitem__(self, key):
"""Access a variable by name."""
return self._data[key]
@property
def data(self):
"""Dict-like view of all computed variable data."""
return self._data
@property
def total_time(self):
if self.set_up_time is None or self.solve_time is None:
return None
return self.set_up_time + self.solve_time
def save(self, filename):
"""Save the entire solution using pickle."""
with open(filename, "wb") as f:
pickle.dump(self, f, pickle.HIGHEST_PROTOCOL)
def get_data_dict(self):
"""Return solution data as a plain ``dict``.
Returns
-------
dict
Mapping of variable names to NumPy arrays.
"""
return dict(self._data)
def save_data(self, filename, to_format="csv", short_names=None):
"""Save solution data to a file.
Parameters
----------
filename : str
The file path to save to.
to_format : str, optional
Output format: ``"csv"`` (default), ``"json"``, ``"pickle"``
or ``"matlab"``.
short_names : dict, optional
Dictionary of shortened names to use when saving.
"""
data = self.get_data_dict()
if short_names:
data = {short_names.get(k, k): v for k, v in data.items()}
if to_format == "csv":
import pandas as _pd
_pd.DataFrame(data).to_csv(filename, index=False)
elif to_format == "json":
with open(filename, "w") as f:
json.dump(data, f, cls=NumpyEncoder)
elif to_format == "pickle":
with open(filename, "wb") as f:
pickle.dump(data, f, pickle.HIGHEST_PROTOCOL)
elif to_format == "matlab":
from scipy.io import savemat as _savemat
_savemat(filename, data)
else:
raise ValueError(
f"Unrecognised format '{to_format}'. "
"Must be 'csv', 'json', 'pickle', or 'matlab'."
)
class EISSolution(SolutionBase):
"""Solution for frequency-domain EIS simulations.
Variables are accessible via dict-style access, e.g.
``solution["Impedance [Ohm]"]``, ``solution["Frequency [Hz]"]``.
Parameters
----------
frequencies : np.ndarray
Frequencies in Hz at which impedance was computed.
impedance : np.ndarray
Complex impedance values at each frequency.
"""
def __init__(self, frequencies, impedance):
super().__init__()
impedance = np.asarray(impedance, dtype=complex)
self._data["Frequency [Hz]"] = np.asarray(frequencies)
self._data["Impedance [Ohm]"] = impedance
self._data["Z_re [Ohm]"] = impedance.real
self._data["Z_im [Ohm]"] = impedance.imag
@property
def frequencies(self):
"""Frequencies in Hz."""
return self._data["Frequency [Hz]"]
@property
def impedance(self):
"""Complex impedance values."""
return self._data["Impedance [Ohm]"]
def get_data_dict(self):
"""Return EIS data as a plain ``dict`` suitable for export.
Returns real-valued columns only (frequencies, real and imaginary
parts of impedance). Use ``solution["Impedance [Ohm]"]`` for the
complex-valued array.
"""
return {
"Frequency [Hz]": self.frequencies,
"Z_re [Ohm]": self._data["Z_re [Ohm]"],
"Z_im [Ohm]": self._data["Z_im [Ohm]"],
}
def nyquist_plot(self, **kwargs):
"""Generate a Nyquist plot from the impedance data.
Parameters
----------
**kwargs
Keyword arguments forwarded to :func:`pybamm.nyquist_plot`.
Returns
-------
fig : matplotlib.figure.Figure or None
ax : matplotlib.axes.Axes
"""
from pybamm.plotting.nyquist_plot import nyquist_plot
return nyquist_plot(self.impedance, **kwargs)
_DEFAULT_SOLUTION_OPTIONS = {
"compile": False,
"cse": True,
"inputs_check": False,
"is_diff_in": [False, True, False],
"is_diff_out": [True],
"regularity_check": False,
"error_on_fail": False,
}
[docs]
class Solution(SolutionBase):
"""
Class containing the solution of, and various attributes associated with, a PyBaMM
model.
Parameters
----------
all_ts : :class:`numpy.array`, size (n,) (or list of these)
A one-dimensional array containing the times at which the solution is evaluated.
A list of times can be provided instead to initialize a solution with
sub-solutions.
all_ys : :class:`numpy.array`, size (m, n) (or list of these)
A two-dimensional array containing the values of the solution. y[i, :] is the
vector of solutions at time t[i].
A list of ys can be provided instead to initialize a solution with
sub-solutions.
all_models : :class:`pybamm.BaseModel`
The model that was used to calculate the solution.
A list of models can be provided instead to initialize a solution with
sub-solutions that have been calculated using those models.
all_inputs : dict (or list of these)
The inputs that were used to calculate the solution
A list of inputs can be provided instead to initialize a solution with
sub-solutions.
t_event : :class:`numpy.array`, size (1,)
A zero-dimensional array containing the time at which the event happens.
y_event : :class:`numpy.array`, size (m,)
A one-dimensional array containing the value of the solution at the time when
the event happens.
termination : str
String to indicate why the solution terminated
all_sensitivities: dict of lists
sensitivities are provided as a dict of {parameter: [sensitivities]} pairs.
all_t_evals : :class:`numpy.array` (or list of these), optional
Evaluation time breakpoints for each segment, one array per entry in
``all_ts``. When provided by the solver these contain the ``t_eval``
times (including merged discontinuity times). If None (the default),
the property returns ``all_ts`` for each segment, since each point
is asssumed to be an evaluation time.
variables_returned: bool
Bool to indicate if `all_ys` contains the full state vector, or is empty because
only requested variables have been returned. True if `output_variables` is used
with a solver, otherwise False.
check_solution: bool
Bool to indicate if the solution should be checked for large y values.
options: dict, optional
Overrides for :meth:`process_casadi_var`. Unspecified keys fall back
to ``_DEFAULT_SOLUTION_OPTIONS`` -- the casadi ``Function`` flags and
the ``compile`` flag (``True`` = AOT, ``False`` = in-process VM).
Solvers should forward the user's selection so observed variables use
the same backend as the integration.
"""
def __init__(
self,
all_ts,
all_ys,
all_models,
all_inputs,
t_event=None,
y_event=None,
termination="final time",
all_sensitivities=None,
all_yps=None,
all_t_evals=None,
variables_returned=False,
check_solution=True,
options=None,
):
if not isinstance(all_ts, list):
all_ts = [all_ts]
self._ensure_sorted_t(all_ts, "all_ts")
if not isinstance(all_ys, list):
all_ys = [all_ys]
if not isinstance(all_models, list):
all_models = [all_models]
self._all_ts = all_ts
self._all_ys = all_ys
self._all_ys_and_sens = all_ys
self._all_models = all_models
if (all_yps is not None) and not isinstance(all_yps, list):
all_yps = [all_yps]
self._all_yps = all_yps
if all_t_evals is None:
all_t_evals = all_ts
else:
if not isinstance(all_t_evals, list):
all_t_evals = [all_t_evals]
self._ensure_t_evals(all_ts=all_ts, all_t_evals=all_t_evals)
self._all_t_evals = all_t_evals
self._user_options = options or {}
self._options = _DEFAULT_SOLUTION_OPTIONS | self._user_options
self.variables_returned = variables_returned
self._observable = self._all_models and all(
model.solution_observable for model in self._all_models
)
# Set up inputs
if not isinstance(all_inputs, list):
all_inputs_copy = dict(all_inputs)
for key, value in all_inputs_copy.items():
if isinstance(value, numbers.Number):
all_inputs_copy[key] = np.array([value])
self.all_inputs = [all_inputs_copy]
else:
self.all_inputs = all_inputs
if all_sensitivities is None:
all_sensitivities = {}
if not isinstance(all_sensitivities, dict):
raise TypeError(
"sensitivities arg needs to be a dict of {parameter: [sensitivities]} "
)
self._all_sensitivities = {}
for key, value in all_sensitivities.items():
if isinstance(value, list):
self._all_sensitivities[key] = value
else:
self._all_sensitivities[key] = [value]
# Events
self._t_event = t_event
self._y_event = y_event
self._termination = termination
self.closest_event_idx = None
super().__init__()
self.integration_time = None
self._all_inputs_stacked = None
self._all_inputs_casadi = None
self._variables = {}
# Add self as sub-solution for compatibility with ProcessedVariable
self._sub_solutions = [self]
# initialize empty cycles
self._cycles = []
# Initialize empty inputs
self._t = None
self._t_eval = None
self._y = None
self._yp = None
self._sensitivities = None
# Initialize empty summary variables
self._summary_variables = None
# Initialise initial start time
self.initial_start_time = None
# Check no ys are too large
if check_solution:
self.check_ys_are_not_too_large()
# Solution now uses CasADi
pybamm.citations.register("Andersson2019")
def has_sensitivities(self) -> bool:
return len(self._all_sensitivities) > 0
@staticmethod
def _ensure_t_evals(all_ts, all_t_evals):
# all_ts is already checked upstream
if all_ts is all_t_evals:
return
if len(all_ts) != len(all_t_evals):
raise ValueError(
"The length of `all_t_evals` must match the length of `all_ts`"
)
Solution._ensure_sorted_t(all_t_evals, "all_t_evals")
for t, t_eval in zip(all_ts, all_t_evals, strict=True):
if len(t) == 0:
continue
if t_eval[0] < t[0] or t_eval[-1] > t[-1]:
raise ValueError(
"The values of the `all_t_evals` segments must be within the same interval as `all_ts`"
)
@staticmethod
def _ensure_sorted_t(ts: list[np.ndarray | casadi.DM | casadi.MX], name: str):
t_f_previous = None
for t in ts:
if len(t) == 0:
continue
if not pybamm.solvers.processed_variable._is_sorted(t):
raise ValueError(
f"Time segment `{name}` must be unique and sorted in increasing order."
)
t_0_current = t[0]
if t_f_previous is not None and t_0_current <= t_f_previous:
raise ValueError(
f"The `{name}` time vector must be strictly increasing across all segments "
"of the sub-solutions."
)
t_f_previous = t[-1]
@property
def t(self) -> np.ndarray:
"""Times at which the solution is evaluated"""
if self._t is None:
self._t = np.concatenate(self.all_ts)
return self._t
@property
def t_eval(self) -> np.ndarray:
"""Evaluation time breakpoints for each segment"""
if self._t_eval is None:
self._t_eval = np.concatenate(self.all_t_evals)
return self._t_eval
@property
def y(self):
"""Values of the solution"""
if self._y is None:
self._y = self._concatenate_states(self.all_ys)
return self._y
@property
def yp(self):
"""Time derivatives of the solution"""
if self.hermite_interpolation and self._yp is None:
self._yp = self._concatenate_states(self.all_yps)
return self._yp
def _concatenate_states(
self, states: list[np.ndarray | casadi.DM | casadi.MX]
) -> np.ndarray | casadi.DM | casadi.MX:
if len(states) == 1:
return states[0]
try:
if isinstance(states[0], casadi.DM | casadi.MX):
return casadi.horzcat(*states)
else:
return np.hstack(states)
except ValueError as error:
raise pybamm.SolverError(
"The solution is made up from different models, so the states cannot be "
"computed explicitly."
) from error
@property
def data(self):
for k, v in self._variables.items():
if k not in self._data:
self._data[k] = v.data
return self._data
@property
def sensitivities(self):
"""Values of the sensitivities. Returns a dict of param_name: np_array"""
if self._sensitivities is None:
self.set_sensitivities()
return self._sensitivities
def set_sensitivities(self):
if not self.has_sensitivities():
self._sensitivities = {}
return
self._sensitivities = {}
for key, sens in self._all_sensitivities.items():
self._sensitivities[key] = np.vstack(sens)
def check_ys_are_not_too_large(self):
# Only check last one so that it doesn't take too long
# We only care about the cases where y is growing too large without any
# restraint, so if y gets large in the middle then comes back down that is ok
t, y, model = self.all_ts[-1], self.all_ys[-1], self.all_models[-1]
t = t[-1]
y = y[:, -1]
max_y_value = pybamm.settings.max_y_value
if y.size == 0 or np.max(y) <= max_y_value:
return
for var in chain(model.rhs.keys(), model.algebraic.keys()):
y_var = self[var.name](t)
if np.any(y_var > max_y_value):
pybamm.logger.error(
f"Solution for '{var}' exceeds the maximum allowed value "
f"of `{max_y_value}. This could be due to "
"incorrect scaling, model formulation, or "
"parameter values. The maximum allowed value is set by "
"'pybammm.settings.max_y_value'."
)
@property
def all_ts(self) -> list[np.ndarray | casadi.DM | casadi.MX]:
return self._all_ts
@property
def all_t_evals(self) -> list[np.ndarray | casadi.DM | casadi.MX]:
return self._all_t_evals
@property
def all_ys(self) -> list[np.ndarray | casadi.DM | casadi.MX]:
return self._all_ys
@property
def all_models(self):
"""Model(s) used for solution"""
return self._all_models
@property
def all_inputs_stacked(self) -> list[np.ndarray]:
if self._all_inputs_stacked is None:
self._all_inputs_stacked = [
np.asarray(list(inp.values())).reshape(-1) for inp in self.all_inputs
]
return self._all_inputs_stacked
@property
def all_inputs_casadi(self) -> list[casadi.DM]:
if self._all_inputs_casadi is None:
self._all_inputs_casadi = [
casadi.vertcat(inp) for inp in self.all_inputs_stacked
]
return self._all_inputs_casadi
@property
def all_yps(self) -> list[np.ndarray | casadi.DM | casadi.MX] | None:
return self._all_yps
@property
def hermite_interpolation(self) -> bool:
return self.all_yps is not None
@property
def t_event(self):
"""Time at which the event happens"""
return self._t_event
@property
def y_event(self):
"""Value of the solution at the time of the event"""
return self._y_event
@property
def termination(self):
"""Reason for termination"""
return self._termination
@termination.setter
def termination(self, value):
"""Updates the reason for termination"""
self._termination = value
@property
def observable(self) -> bool:
return self._observable
def _check_observable(self):
if self._observable:
return
# Collect unique reasons from all models
unique_reasons = set(
model.solution_observable_status
for model in self.all_models
if not model.solution_observable
)
if unique_reasons:
reasons_str = ", ".join(sorted(unique_reasons))
raise ValueError(f"Solution is not observable: {reasons_str}")
@property
def user_options(self) -> dict:
return self._user_options
@property
def options(self) -> dict:
return self._options
[docs]
@cached_property
def first_state(self):
"""
A Solution object that only contains the first state. This is faster to evaluate
than the full solution when only the first state is needed (e.g. to initialize
a model with the solution)
"""
sensitivities = {}
n_states = self.all_models[0].len_rhs_and_alg
for key in self._all_sensitivities:
sensitivities[key] = self._all_sensitivities[key][0][-n_states:, :]
if self.all_yps is None:
all_yps = None
else:
all_yps = self.all_yps[0][:, :1]
if not self.variables_returned:
all_ys = self.all_ys[0][:, :1]
else:
# Get first state from initial conditions as all_ys is empty
all_ys = self.all_models[0].y0full[0].reshape(-1, 1)
new_sol = Solution(
self.all_ts[0][:1],
all_ys,
self.all_models[:1],
self.all_inputs[:1],
None,
None,
"final time",
all_sensitivities=sensitivities,
all_yps=all_yps,
options=self.user_options,
)
new_sol._all_inputs_stacked = self.all_inputs_stacked[:1]
new_sol._all_inputs_casadi = self.all_inputs_casadi[:1]
new_sol._sub_solutions = self.sub_solutions[:1]
new_sol.solve_time = 0
new_sol.integration_time = 0
new_sol.set_up_time = 0
return new_sol
[docs]
@cached_property
def last_state(self):
"""
A Solution object that only contains the final state. This is faster to evaluate
than the full solution when only the final state is needed (e.g. to initialize
a model with the solution)
"""
sensitivities = {}
n_states = self.all_models[-1].len_rhs_and_alg
for key in self._all_sensitivities:
sensitivities[key] = self._all_sensitivities[key][-1][-n_states:, :]
if self.all_yps is None:
all_yps = None
else:
all_yps = self.all_yps[-1][:, -1:]
if not self.variables_returned:
all_ys = self.all_ys[-1][:, -1:]
else:
# Get last state from y_event as all_ys is empty
all_ys = self.y_event.reshape(len(self.y_event), 1)
new_sol = Solution(
self.all_ts[-1][-1:],
all_ys,
self.all_models[-1:],
self.all_inputs[-1:],
self.t_event,
self.y_event,
self.termination,
all_sensitivities=sensitivities,
all_yps=all_yps,
options=self.user_options,
)
new_sol._all_inputs_stacked = self.all_inputs_stacked[-1:]
new_sol._all_inputs_casadi = self.all_inputs_casadi[-1:]
new_sol._sub_solutions = self.sub_solutions[-1:]
new_sol.solve_time = 0
new_sol.integration_time = 0
new_sol.set_up_time = 0
return new_sol
@property
def cycles(self):
return self._cycles
@cycles.setter
def cycles(self, cycles):
self._cycles = cycles
@property
def summary_variables(self):
return self._summary_variables
@property
def initial_start_time(self):
return self._initial_start_time
@initial_start_time.setter
def initial_start_time(self, value):
"""Updates the initial start time of the experiment"""
self._initial_start_time = value
def update_summary_variables(self, all_summary_variables):
self.all_summary_variables = all_summary_variables
self._summary_variables = pybamm.SummaryVariables(
self, cycle_summary_variables=all_summary_variables
)
[docs]
def update(self, variables: str | list[str]):
"""Add ProcessedVariables to the dictionary of variables in the solution"""
# Single variable
if isinstance(variables, str):
variables = [variables]
# Process
for variable in variables:
self._update_variable(variable)
def _update_model_variable(
self,
model: pybamm.BaseModel,
var_pybamm: pybamm.Symbol,
time_integral: pybamm.ProcessedVariableTimeIntegral | None,
inputs: dict,
ys_shape: tuple,
cache_key,
):
_var_casadi = model._variables_casadi.get(cache_key)
if _var_casadi is not None:
return _var_casadi, var_pybamm, time_integral
var_casadi, var_pybamm, time_integral = self._convert_to_casadi(
var_pybamm, inputs, ys_shape
)
# Only cache if it's not a time integral
if time_integral is None:
model._variables_casadi[cache_key] = var_casadi
return var_casadi, var_pybamm, time_integral
def _update_variable(self, name: str):
time_integral = None
pybamm.logger.debug(f"Post-processing {name}")
# Iterate through all models, some may be in the list several times and
# therefore only get set up once
vars_pybamm = [
model.get_processed_variable_or_event(name) for model in self.all_models
]
vars_casadi = [None] * len(self.all_models)
for i, (model, ys, inputs) in enumerate(
zip(self.all_models, self.all_ys, self.all_inputs, strict=True)
):
_var_pybamm = vars_pybamm[i]
if self.variables_returned and _var_pybamm.has_symbol_of_classes(
pybamm.expression_tree.state_vector.StateVector
):
raise KeyError(
f"Cannot process variable '{name}' as it was not part of the "
"solve. Please re-run the solve with `output_variables` set to "
"include this variable."
)
var_casadi, var_pybamm, time_integral = self._update_model_variable(
model,
_var_pybamm,
inputs=inputs,
ys_shape=ys.shape,
time_integral=time_integral,
cache_key=name,
)
vars_pybamm[i] = var_pybamm
vars_casadi[i] = var_casadi
var = pybamm.process_variable(
name, vars_pybamm, vars_casadi, self, time_integral=time_integral
)
self._variables[name] = var
[docs]
def observe(self, symbol: pybamm.Symbol) -> pybamm.ProcessedVariable:
"""
Observe a `pybamm.Symbol` object from the solution.
Parameters
----------
symbol : pybamm.Symbol
The symbol to observe.
Returns
-------
:class:`pybamm.ProcessedVariable`
The observed variable.
"""
self._check_observable()
symbol = pybamm.convert_to_symbol(symbol)
name = str(symbol.id)
value = self._variables.get(name, None)
if value is not None:
return value
for model in self.all_models:
model.process_and_register_variable(name=name, symbol=symbol)
return self[name]
def _convert_to_casadi(self, var_pybamm, inputs, ys_shape):
time_integral = pybamm.ProcessedVariableTimeIntegral.from_pybamm_var(
var_pybamm, ys_shape[0]
)
if time_integral is not None:
var_pybamm = time_integral.sum_node.child
if time_integral.post_sum_node is not None:
time_integral.post_sum = self.process_casadi_var(
time_integral.post_sum_node,
inputs,
ys_shape,
)
var_casadi = self.process_casadi_var(
var_pybamm,
inputs,
ys_shape,
)
return var_casadi, var_pybamm, time_integral
_CASADI_FUNCTION_OPTIONS = [
"cse",
"inputs_check",
"is_diff_in",
"is_diff_out",
"regularity_check",
"error_on_fail",
]
def process_casadi_var(self, var_pybamm, inputs, ys_shape):
t_MX = casadi.MX.sym("t")
y_MX = casadi.MX.sym("y", ys_shape[0])
total_input_size = sum(v.size for v in inputs.values())
inputs_MX = casadi.MX.sym("p", total_input_size)
inputs_MX_dict = {}
offset = 0
for key, value in inputs.items():
n = value.size
inputs_MX_dict[key] = inputs_MX[offset : offset + n]
offset += n
var_sym = var_pybamm.to_casadi(t_MX, y_MX, inputs=inputs_MX_dict)
opts = {opt: self.options[opt] for opt in self._CASADI_FUNCTION_OPTIONS}
# Casadi has a bug where it does not correctly handle arrays with
# zeros padded at the beginning or end. To avoid this, we add and
# subtract the same number to the variable to reinforce the
# variable bounds. This does not affect the answer
epsilon = 1.0
var_sym = (var_sym - epsilon) + epsilon
var_casadi = casadi.Function(
"variable",
[t_MX, y_MX, inputs_MX],
[var_sym],
opts,
)
# Skip the SX expansion on the aot path: it only exists to give the
# vm interpreter a faster bytecode form, and compiled code doesn't
# need it.
if self._options["compile"]:
return aot_compile(var_casadi)
# Some variables, like interpolants, cannot be expanded
try:
var_casadi_out = var_casadi.expand()
except RuntimeError as error:
if "'eval_sx' not defined for" not in str(error):
raise error # pragma: no cover
var_casadi_out = var_casadi
return var_casadi_out
def __getitem__(self, key):
"""Read a variable from the solution. Variables are created 'just in time', i.e.
only when they are called.
Parameters
----------
key : str
The name of the variable
Returns
-------
:class:`pybamm.ProcessedVariable` or :class:`pybamm.ProcessedVariableComputed`
A variable that can be evaluated at any time or spatial point. The
underlying data for this variable is available in its attribute ".data"
"""
value = self._variables.get(key)
if value is not None:
return value
# otherwise create it, save it and then return it
self.update(key)
return self._variables[key]
[docs]
def plot(self, output_variables=None, **kwargs):
"""
A method to quickly plot the outputs of the solution. Creates a
:class:`pybamm.QuickPlot` object (with keyword arguments 'kwargs') and
then calls :meth:`pybamm.QuickPlot.dynamic_plot`.
Parameters
----------
output_variables: list, optional
A list of the variables to plot.
**kwargs
Additional keyword arguments passed to
:meth:`pybamm.QuickPlot.dynamic_plot`.
For a list of all possible keyword arguments see :class:`pybamm.QuickPlot`.
"""
return pybamm.dynamic_plot(self, output_variables=output_variables, **kwargs)
[docs]
def get_data_dict(self, variables=None, short_names=None, cycles_and_steps=True):
"""
Construct a (standard python) dictionary of the solution data containing the
variables in `variables`. If `variables` is None then all variables are
returned. Any variable names in short_names are replaced with the corresponding
short name.
If the solution has cycles, then the cycle numbers and step numbers are also
returned in the dictionary.
Parameters
----------
variables : list, optional
List of variables to return. If None, returns all variables in solution.data
short_names : dict, optional
Dictionary of shortened names to use when saving.
cycles_and_steps : bool, optional
Whether to include the cycle numbers and step numbers in the dictionary
Returns
-------
dict
A dictionary of the solution data
"""
if variables is None:
# variables not explicitly provided -> save all variables that have been
# computed
data_long_names = self.data
else:
if isinstance(variables, str):
variables = [variables]
# otherwise, save only the variables specified
data_long_names = {}
for name in variables:
data_long_names[name] = self[name].data
if len(data_long_names) == 0:
raise ValueError(
"""
Solution does not have any data. Please provide a list of variables
to save.
"""
)
# Use any short names if provided
data_short_names = {}
short_names = short_names or {}
for name, var in data_long_names.items():
name = short_names.get(name, name) # return name if no short name
data_short_names[name] = var
# Save cycle number and step number if the solution has them
if cycles_and_steps and len(self.cycles) > 0:
data_short_names["Cycle"] = np.array([])
data_short_names["Step"] = np.array([])
for i, cycle in enumerate(self.cycles):
data_short_names["Cycle"] = np.concatenate(
[data_short_names["Cycle"], i * np.ones_like(cycle.t)]
)
for j, step in enumerate(cycle.steps):
data_short_names["Step"] = np.concatenate(
[data_short_names["Step"], j * np.ones_like(step.t)]
)
return data_short_names
[docs]
def save_data(
self, filename=None, variables=None, to_format="pickle", short_names=None
):
"""
Save solution data only (raw arrays)
Parameters
----------
filename : str, optional
The name of the file to save data to. If None, then a str is returned
variables : list, optional
List of variables to save. If None, saves all of the variables that have
been created so far
to_format : str, optional
The format to save to. Options are:
- 'pickle' (default): creates a pickle file with the data dictionary
- 'matlab': creates a .mat file, for loading in matlab
- 'csv': creates a csv file (0D variables only)
- 'json': creates a json file
short_names : dict, optional
Dictionary of shortened names to use when saving. This may be necessary when
saving to MATLAB, since no spaces or special characters are allowed in
MATLAB variable names. Note that not all the variables need to be given
a short name.
Returns
-------
data : str, optional
str if 'csv' or 'json' is chosen and filename is None, otherwise None
"""
data = self.get_data_dict(variables=variables, short_names=short_names)
if to_format == "pickle":
if filename is None:
raise ValueError("pickle format must be written to a file")
with open(filename, "wb") as f:
pickle.dump(data, f, pickle.HIGHEST_PROTOCOL)
elif to_format == "matlab":
if filename is None:
raise ValueError("matlab format must be written to a file")
# Check all the variable names only contain a-z, A-Z or _ or numbers
for name in data.keys():
# Check the string only contains the following ASCII:
# a-z (97-122)
# A-Z (65-90)
# _ (95)
# 0-9 (48-57) but not in the first position
for i, s in enumerate(name):
if not (
97 <= ord(s) <= 122
or 65 <= ord(s) <= 90
or ord(s) == 95
or (i > 0 and 48 <= ord(s) <= 57)
):
raise ValueError(
f"Invalid character '{s}' found in '{name}'. "
+ "MATLAB variable names must only contain a-z, A-Z, _, "
"or 0-9 (except the first position). "
"Use the 'short_names' argument to pass an alternative "
"variable name, e.g. \n\n"
"\tsolution.save_data(filename, "
"['Electrolyte concentration'], to_format='matlab, "
"short_names={'Electrolyte concentration': 'c_e'})"
)
savemat(filename, data)
elif to_format == "csv":
for name, var in data.items():
if var.ndim >= 2:
raise ValueError(
f"only 0D variables can be saved to csv, but '{name}' is {var.ndim - 1}D"
)
df = pd.DataFrame(data)
return df.to_csv(filename, index=False)
elif to_format == "json":
if filename is None:
return json.dumps(data, cls=NumpyEncoder)
else:
with open(filename, "w") as outfile:
json.dump(data, outfile, cls=NumpyEncoder)
else:
raise ValueError(f"format '{to_format}' not recognised")
@property
def sub_solutions(self):
"""List of sub solutions that have been
concatenated to form the full solution"""
return self._sub_solutions
def __add__(self, other):
"""Adds two solutions together, e.g. when stepping"""
if other is None or isinstance(other, EmptySolution):
return self.copy()
if not isinstance(other, Solution):
raise pybamm.SolverError(
"Only a Solution or None can be added to a Solution"
)
# Special case: new solution only has one timestep and it is already in the
# existing solution. In this case, return a copy of the existing solution
if (
len(other.all_ts) == 1
and len(other.all_ts[0]) == 1
and other.all_ts[0][0] == self.all_ts[-1][-1]
):
new_sol = self.copy()
# Update termination using the latter solution
new_sol._termination = other.termination
new_sol._t_event = other._t_event
new_sol._y_event = other._y_event
return new_sol
# Update list of sub-solutions
hermite_interpolation = (
other.hermite_interpolation and self.hermite_interpolation
)
if other.all_ts[0][0] == self.all_ts[-1][-1]:
# Skip first time step if it is repeated
all_ts = [*self.all_ts, other.all_ts[0][1:], *other.all_ts[1:]]
all_ys = [*self.all_ys, other.all_ys[0][:, 1:], *other.all_ys[1:]]
if hermite_interpolation:
all_yps = [*self.all_yps, other.all_yps[0][:, 1:], *other.all_yps[1:]]
if self._all_t_evals is not None and other._all_t_evals is not None:
all_t_evals = [
*self._all_t_evals,
other._all_t_evals[0][1:],
*other._all_t_evals[1:],
]
else:
all_t_evals = None
else:
all_ts = self.all_ts + other.all_ts
all_ys = self.all_ys + other.all_ys
if hermite_interpolation:
all_yps = self.all_yps + other.all_yps
if self._all_t_evals is not None and other._all_t_evals is not None:
all_t_evals = self._all_t_evals + other._all_t_evals
else:
all_t_evals = None
if not hermite_interpolation:
all_yps = None
# sensitivities are a dict of {parameter: [sensitivities]}
# we can assume that the keys are the same for both solutions
all_sensitivities = self._all_sensitivities
for key in other._all_sensitivities:
all_sensitivities[key] = (
all_sensitivities[key] + other._all_sensitivities[key]
)
options = self.user_options | other.user_options
new_sol = Solution(
all_ts,
all_ys,
self.all_models + other.all_models,
self.all_inputs + other.all_inputs,
other.t_event,
other.y_event,
other.termination,
all_sensitivities=all_sensitivities,
all_yps=all_yps,
all_t_evals=all_t_evals,
variables_returned=other.variables_returned,
options=options,
)
new_sol.closest_event_idx = other.closest_event_idx
new_sol._all_inputs_stacked = self.all_inputs_stacked + other.all_inputs_stacked
new_sol._all_inputs_casadi = self.all_inputs_casadi + other.all_inputs_casadi
# Add timers (if available)
for attr in ["solve_time", "integration_time", "set_up_time"]:
if (
getattr(self, attr, None) is not None
and getattr(other, attr, None) is not None
):
setattr(new_sol, attr, getattr(self, attr) + getattr(other, attr))
# Set sub_solutions
new_sol._sub_solutions = self.sub_solutions + other.sub_solutions
# update variables which were derived at the solver stage
if any([self.variables_returned, other.variables_returned]):
vars = {*self._variables.keys(), *other._variables.keys()}
new_sol._variables = {v: self[v].update(other[v], new_sol) for v in vars}
new_sol.variables_returned = True
return new_sol
def __radd__(self, other):
return self.__add__(other)
def copy(self):
new_sol = self.__class__(
self.all_ts,
self.all_ys,
self.all_models,
self.all_inputs,
self.t_event,
self.y_event,
self.termination,
all_sensitivities=self._all_sensitivities,
all_yps=self.all_yps,
all_t_evals=self._all_t_evals,
variables_returned=self.variables_returned,
options=self.user_options,
)
new_sol._all_inputs_stacked = self.all_inputs_stacked
new_sol._all_inputs_casadi = self.all_inputs_casadi
new_sol._sub_solutions = self.sub_solutions
new_sol.closest_event_idx = self.closest_event_idx
new_sol.solve_time = self.solve_time
new_sol.integration_time = self.integration_time
new_sol.set_up_time = self.set_up_time
# copy over variables which were derived at the solver stage
if self._variables and all(
isinstance(v, pybamm.ProcessedVariableComputed)
for v in self._variables.values()
):
new_sol._variables = self._variables.copy()
return new_sol
[docs]
def plot_voltage_components(
self,
ax=None,
show_legend=True,
split_by_electrode=False,
show_plot=True,
**kwargs_fill,
):
"""
Generate a plot showing the component overpotentials that make up the voltage
Parameters
----------
ax : matplotlib Axis, optional
The axis on which to put the plot. If None, a new figure and axis is created.
show_legend : bool, optional
Whether to display the legend. Default is True.
split_by_electrode : bool, optional
Whether to show the overpotentials for the negative and positive electrodes
separately. Default is False.
show_plot : bool, optional
Whether to show the plots. Default is True. Set to False if you want to
only display the plot after plt.show() has been called.
kwargs_fill
Keyword arguments, passed to ax.fill_between.
"""
return pybamm.plot_voltage_components(
self,
ax=ax,
show_legend=show_legend,
split_by_electrode=split_by_electrode,
show_plot=show_plot,
**kwargs_fill,
)
class EmptySolution:
def __init__(self, termination=None, t=None):
self.termination = termination
if t is None:
t = np.array([0])
elif isinstance(t, numbers.Number):
t = np.array([t])
self.t = t
def __add__(self, other):
if isinstance(other, EmptySolution | Solution):
return other.copy()
def __radd__(self, other):
if other is None:
return self.copy()
def copy(self):
return EmptySolution(termination=self.termination, t=self.t)
def make_cycle_solution(
step_solutions, esoh_solver=None, save_this_cycle=True, inputs=None
):
"""
Function to create a Solution for an entire cycle, and associated summary variables
Parameters
----------
step_solutions : list of :class:`Solution`
Step solutions that form the entire cycle
esoh_solver : :class:`pybamm.lithium_ion.ElectrodeSOHSolver`
Solver to calculate electrode SOH (eSOH) variables. If `None` (default)
then only summary variables that do not require the eSOH calculation
are calculated. See :footcite:t:`Mohtat2019` for more details on eSOH variables.
save_this_cycle : bool, optional
Whether to save the entire cycle variables or just the summary variables.
Default True
inputs : dict
User inputs for the summary variables.
Returns
-------
cycle_solution : :class:`pybamm.Solution` or None
The Solution object for this cycle, or None (if save_this_cycle is False)
cycle_summary_variables : dict
Dictionary of summary variables for this cycle
cycle_first_state : Solution
First state of the cycle.
"""
sum_sols = step_solutions[0].copy()
for step_solution in step_solutions[1:]:
sum_sols = sum_sols + step_solution
cycle_solution = Solution(
sum_sols.all_ts,
sum_sols.all_ys,
sum_sols.all_models,
sum_sols.all_inputs,
sum_sols.t_event,
sum_sols.y_event,
sum_sols.termination,
all_sensitivities=sum_sols._all_sensitivities,
all_yps=sum_sols.all_yps,
all_t_evals=sum_sols._all_t_evals,
variables_returned=sum_sols.variables_returned,
options=sum_sols.user_options,
)
if sum_sols.variables_returned:
cycle_solution._variables = sum_sols._variables
cycle_solution._all_inputs_stacked = sum_sols.all_inputs_stacked
cycle_solution._all_inputs_casadi = sum_sols.all_inputs_casadi
cycle_solution._sub_solutions = sum_sols.sub_solutions
cycle_solution.solve_time = sum_sols.solve_time
cycle_solution.integration_time = sum_sols.integration_time
cycle_solution.set_up_time = sum_sols.set_up_time
cycle_solution.steps = step_solutions
cycle_summary_variables = pybamm.SummaryVariables(
cycle_solution, esoh_solver=esoh_solver, user_inputs=inputs
)
cycle_first_state = cycle_solution.first_state
if save_this_cycle:
cycle_solution.cycle_summary_variables = cycle_summary_variables
else:
cycle_solution = None
return cycle_solution, cycle_summary_variables, cycle_first_state
