Fix indexing bug with infeasible experiments for IDAKLUSolver (#4541)

* fix interp indexing

Co-Authored-By: Pip Liggins <[email protected]>

* simplify indexing

---------

Co-authored-by: Pip Liggins <[email protected]>

Author: MarcBerliner

CHANGELOG.md

@@ -20,7 +20,7 @@
 - Removed the `start_step_offset` setting and disabled minimum `dt` warnings for drive cycles with the (`IDAKLUSolver`). ([#4416](https://github.com/pybamm-team/PyBaMM/pull/4416))
 
 ## Bug Fixes
-
+- Fixed bug in post-processing solutions with infeasible experiments using the (`IDAKLUSolver`). ([#4541](https://github.com/pybamm-team/PyBaMM/pull/4541))
 - Disabled IREE on MacOS due to compatibility issues and added the CasADI
   path to the environment to resolve issues on MacOS and Linux. Windows
   users may still experience issues with interpolation. ([#4528](https://github.com/pybamm-team/PyBaMM/pull/4528))

src/pybamm/solvers/c_solvers/idaklu/observe.cpp

@@ -100,6 +100,11 @@ class TimeSeriesInterpolator {
         ) {
         for (size_t i = 0; i < ts_data_np.size(); i++) {
             const auto& t_data = ts_data_np[i].unchecked<1>();
+            // Continue if there is no data
+            if (t_data.size() == 0) {
+                continue;
+            }
+
             const realtype t_data_final = t_data(t_data.size() - 1);
             realtype t_interp_next = t_interp(i_interp);
             // Continue if the next interpolation point is beyond the final data point
@@ -227,6 +232,10 @@ class TimeSeriesProcessor {
         int i_entries = 0;
         for (size_t i = 0; i < ts.size(); i++) {
             const auto& t = ts[i].unchecked<1>();
+            // Continue if there is no data
+            if (t.size() == 0) {
+                continue;
+            }
             const auto& y = ys[i].unchecked<2>();
             const auto input = inputs[i].data();
             const auto func = *funcs[i];

src/pybamm/solvers/processed_variable.py

@@ -133,19 +133,22 @@ def _setup_cpp_inputs(self, t, full_range):
         ys = self.all_ys
         yps = self.all_yps
         inputs = self.all_inputs_casadi
-        # Find the indices of the time points to observe
-        if full_range:
-            idxs = range(len(ts))
-        else:
-            idxs = _find_ts_indices(ts, t)
 
-        if isinstance(idxs, list):
-            # Extract the time points and inputs
-            ts = [ts[idx] for idx in idxs]
-            ys = [ys[idx] for idx in idxs]
-            if self.hermite_interpolation:
-                yps = [yps[idx] for idx in idxs]
-            inputs = [self.all_inputs_casadi[idx] for idx in idxs]
+        # Remove all empty ts
+        idxs = np.where([ti.size > 0 for ti in ts])[0]
+
+        # Find the indices of the time points to observe
+        if not full_range:
+            ts_nonempty = [ts[idx] for idx in idxs]
+            idxs_subset = _find_ts_indices(ts_nonempty, t)
+            idxs = idxs[idxs_subset]
+
+        # Extract the time points and inputs
+        ts = [ts[idx] for idx in idxs]
+        ys = [ys[idx] for idx in idxs]
+        if self.hermite_interpolation:
+            yps = [yps[idx] for idx in idxs]
+        inputs = [self.all_inputs_casadi[idx] for idx in idxs]
 
         is_f_contiguous = _is_f_contiguous(ys)
 
@@ -977,8 +980,4 @@ def _find_ts_indices(ts, t):
     if (t[-1] > ts[-1][-1]) and (len(indices) == 0 or indices[-1] != len(ts) - 1):
         indices.append(len(ts) - 1)
 
-    if len(indices) == len(ts):
-        # All indices are included
-        return range(len(ts))
-
     return indices

src/pybamm/solvers/solution.py

@@ -580,11 +580,15 @@ def _update_variable(self, variable):
         # Iterate through all models, some may be in the list several times and
         # therefore only get set up once
         vars_casadi = []
-        for i, (model, ys, inputs, var_pybamm) in enumerate(
-            zip(self.all_models, self.all_ys, self.all_inputs, vars_pybamm)
+        for i, (model, ts, ys, inputs, var_pybamm) in enumerate(
+            zip(self.all_models, self.all_ts, self.all_ys, self.all_inputs, vars_pybamm)
         ):
-            if ys.size == 0 and var_pybamm.has_symbol_of_classes(
-                pybamm.expression_tree.state_vector.StateVector
+            if (
+                ys.size == 0
+                and var_pybamm.has_symbol_of_classes(
+                    pybamm.expression_tree.state_vector.StateVector
+                )
+                and not ts.size == 0
             ):
                 raise KeyError(
                     f"Cannot process variable '{variable}' as it was not part of the "