Semiquant: Add spline knots to the optimization result

doc/example/semiquantitative_data.ipynb

@@ -431,7 +431,10 @@
     "    )\n",
     "\n",
     "    plot_splines_from_inner_result(\n",
-    "        inner_problem, inner_solvers[minimal_diff], results[minimal_diff]\n",
+    "        inner_problem,\n",
+    "        inner_solvers[minimal_diff],\n",
+    "        results[minimal_diff],\n",
+    "        sim=[simulation],\n",
     "    )\n",
     "    plt.show()"
    ]
@@ -467,7 +470,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.2"
+   "version": "3.10.10"
   },
   "vscode": {
    "interpreter": {

pypesto/C.py

@@ -89,7 +89,6 @@ class EnsembleType(Enum):
 INNER_PARAMETERS = "inner_parameters"
 INNER_RDATAS = "inner_rdatas"
 PARAMETER_TYPE = "parameterType"
-X_INNER_OPT = "x_inner_opt"
 RELATIVE = "relative"
 
 
@@ -207,6 +206,7 @@ class InnerParameterType(str, Enum):
 MIN_SIM_RANGE = 1e-16
 
 SPLINE_PAR_TYPE = "spline"
+SPLINE_KNOTS = "spline_knots"
 N_SPLINE_PARS = "n_spline_pars"
 DATAPOINTS = "datapoints"
 MIN_DATAPOINT = "min_datapoint"

pypesto/hierarchical/inner_calculator_collector.py

@@ -31,9 +31,10 @@
     RES,
     SEMIQUANTITATIVE,
     SPLINE_APPROXIMATION_OPTIONS,
+    SPLINE_KNOTS,
     SPLINE_RATIO,
     SRES,
-    X_INNER_OPT,
+    InnerParameterType,
     ModeType,
 )
 from ..objective.amici.amici_calculator import AmiciCalculator
@@ -107,6 +108,7 @@ def __init__(
         self.quantitative_data_mask = self._get_quantitative_data_mask(edatas)
 
         self._known_least_squares_safe = False
+        self.semiquant_observable_ids = None
 
     def initialize(self):
         """Initialize."""
@@ -177,6 +179,12 @@ def construct_inner_calculators(
                 semiquant_problem.get_noise_dummy_values(scaled=True)
             )
             self.inner_calculators.append(semiquant_calculator)
+            self.semiquant_observable_ids = [
+                model.getObservableIds()[group - 1]
+                for group in semiquant_problem.get_groups_for_xs(
+                    InnerParameterType.SPLINE
+                )
+            ]
 
         if self.data_types - {
             RELATIVE,
@@ -382,7 +390,7 @@ def __call__(
         nllh, snllh, s2nllh, chi2, res, sres = init_return_values(
             sensi_orders, mode, dim
         )
-        all_inner_pars = {}
+        spline_knots = None
         interpretable_inner_pars = []
 
         # set order in solver
@@ -421,7 +429,7 @@ def __call__(
                 RES: res,
                 SRES: sres,
                 RDATAS: rdatas,
-                X_INNER_OPT: all_inner_pars,
+                SPLINE_KNOTS: None,
                 INNER_PARAMETERS: None,
             }
             ret[FVAL] = np.inf
@@ -473,10 +481,11 @@ def __call__(
             if 1 in sensi_orders:
                 snllh += inner_result[GRAD]
 
-            all_inner_pars.update(inner_result[X_INNER_OPT])
             inner_pars = inner_result.get(INNER_PARAMETERS)
             if inner_pars is not None:
                 interpretable_inner_pars.extend(inner_pars)
+            if SPLINE_KNOTS in inner_result:
+                spline_knots = inner_result[SPLINE_KNOTS]
 
         # add the quantitative data contribution
         if self.quantitative_data_mask is not None:
@@ -507,7 +516,7 @@ def __call__(
         # Add inner parameters to return dict
         # only if the objective value improved.
         if ret[FVAL] < self.best_fval:
-            ret[X_INNER_OPT] = all_inner_pars
+            ret[SPLINE_KNOTS] = spline_knots
             ret[INNER_PARAMETERS] = (
                 interpretable_inner_pars
                 if len(interpretable_inner_pars) > 0

pypesto/hierarchical/ordinal/calculator.py

@@ -17,7 +17,6 @@
     RDATAS,
     RES,
     SRES,
-    X_INNER_OPT,
 )
 from ...objective.amici.amici_calculator import (
     AmiciCalculator,
@@ -126,7 +125,7 @@ def __call__(
         Returns
         -------
         inner_result:
-            A dict containing the calculation results: FVAL, GRAD, RDATAS and X_INNER_OPT.
+            A dict containing the calculation results: FVAL, GRAD, RDATAS.
         """
         if mode == MODE_RES:
             raise ValueError(
@@ -178,7 +177,6 @@ def __call__(
             RES: res,
             SRES: sres,
             RDATAS: rdatas,
-            X_INNER_OPT: self.inner_problem.get_inner_parameter_dictionary(),
         }
 
         # if any amici simulation failed, it's unlikely we can compute
@@ -201,13 +199,9 @@ def __call__(
         inner_result[FVAL] = self.inner_solver.calculate_obj_function(
             x_inner_opt
         )
-        inner_result[
-            X_INNER_OPT
-        ] = self.inner_problem.get_inner_parameter_dictionary()
 
         # calculate analytical gradients if requested
         if sensi_order > 0:
-            # print([opt['fun'] for opt in x_inner_opt])
             sy = [rdata[AMICI_SY] for rdata in rdatas]
             ssigma = [rdata[AMICI_SSIGMAY] for rdata in rdatas]
             inner_result[GRAD] = self.inner_solver.calculate_gradients(

pypesto/hierarchical/relative/calculator.py

@@ -25,7 +25,6 @@
     RDATAS,
     RES,
     SRES,
-    X_INNER_OPT,
     ModeType,
 )
 from ...objective.amici.amici_calculator import (
@@ -123,7 +122,7 @@ def __call__(
         Returns
         -------
         inner_result:
-            A dict containing the calculation results: FVAL, GRAD, RDATAS and X_INNER_OPT.
+            A dict containing the calculation results: FVAL, GRAD, RDATAS and INNER_PARAMETERS.
         """
         if not self.inner_problem.check_edatas(edatas=edatas):
             raise ValueError(
@@ -164,8 +163,6 @@ def __call__(
                 rdatas=rdatas,
             )
 
-        inner_result[X_INNER_OPT] = {}
-
         inner_result[INNER_PARAMETERS] = (
             np.array(
                 [

pypesto/hierarchical/semiquantitative/calculator.py

@@ -15,8 +15,8 @@
     MODE_RES,
     RDATAS,
     RES,
+    SPLINE_KNOTS,
     SRES,
-    X_INNER_OPT,
 )
 from ...objective.amici.amici_calculator import (
     AmiciCalculator,
@@ -119,7 +119,8 @@ def __call__(
         Returns
         -------
         inner_result:
-            A dict containing the calculation results: FVAL, GRAD, RDATAS and X_INNER_OPT.
+            A dict containing the calculation results: FVAL, GRAD, RDATAS,
+            INNER_PARAMETERS, and SPLINE_KNOTS.
         """
         if mode == MODE_RES:
             raise ValueError(
@@ -175,7 +176,6 @@ def __call__(
             RES: res,
             SRES: sres,
             RDATAS: rdatas,
-            X_INNER_OPT: self.inner_problem.get_inner_parameter_dictionary(),
         }
 
         # if any amici simulation failed, it's unlikely we can compute
@@ -198,9 +198,7 @@ def __call__(
         inner_result[FVAL] = self.inner_solver.calculate_obj_function(
             x_inner_opt
         )
-        inner_result[
-            X_INNER_OPT
-        ] = self.inner_problem.get_inner_parameter_dictionary()
+        inner_result[SPLINE_KNOTS] = self.inner_problem.get_spline_knots()
 
         inner_result[
             INNER_PARAMETERS

pypesto/hierarchical/semiquantitative/problem.py

@@ -205,6 +205,50 @@ def get_inner_parameter_dictionary(self) -> dict:
             inner_par_dict[x_id] = x.value
         return inner_par_dict
 
+    def get_spline_knots(
+        self,
+    ) -> list[list[np.ndarray[float], np.ndarray[float]]]:
+        """Get spline knots of all semiquantitative observables.
+
+        Returns
+        -------
+        list[list[np.ndarray[float], np.ndarray[float]]]
+            A list of lists with two arrays. Each list in the first level corresponds
+            to a semiquantitative observable. Each of these lists contains two arrays:
+            the first array contains the spline bases, the second array contains the
+            spline knot values. The ordering of the observable lists is the same
+            as in `pypesto.problem.hierarchical.semiquant_observable_ids`.
+        """
+        # We need the solver only for the rescaling function.
+        from .solver import SemiquantInnerSolver
+
+        all_spline_knots = []
+
+        for group in self.get_groups_for_xs(InnerParameterType.SPLINE):
+            group_dict = self.groups[group]
+            n_spline_pars = group_dict[N_SPLINE_PARS]
+            n_data_points = group_dict[NUM_DATAPOINTS]
+
+            inner_pars = np.array(
+                [x.value for x in self.get_xs_for_group(group)]
+            )
+
+            # Utility matrix for the spline knot calculation
+            lower_trian = np.tril(np.ones((n_spline_pars, n_spline_pars)))
+            knot_values = np.dot(lower_trian, inner_pars)
+
+            _, knot_bases, _ = SemiquantInnerSolver._rescale_spline_bases(
+                sim_all=group_dict[CURRENT_SIMULATION],
+                N=n_spline_pars,
+                K=n_data_points,
+            )
+
+            spline_knots_for_observable = [knot_bases, knot_values]
+
+            all_spline_knots.append(spline_knots_for_observable)
+
+        return all_spline_knots
+
     def get_measurements_for_group(self, gr) -> np.ndarray:
         """Get measurements for a group."""
         # Taking the ixs of first inner parameter since

pypesto/hierarchical/semiquantitative/solver.py

@@ -440,7 +440,8 @@ def inner_gradient_wrapper(x):
 
         return results
 
-    def _rescale_spline_bases(self, sim_all: np.ndarray, N: int, K: int):
+    @staticmethod
+    def _rescale_spline_bases(sim_all: np.ndarray, N: int, K: int):
         """Rescale the spline bases.
 
         Before the optimization of the spline parameters, we have to fix the

pypesto/history/csv.py

@@ -20,7 +20,6 @@
     RES,
     SRES,
     TIME,
-    X_INNER_OPT,
     ModeType,
     X,
 )
@@ -155,10 +154,6 @@ def _update_trace(
             else:
                 row[(var, np.nan)] = np.nan
 
-        if X_INNER_OPT in result:
-            for x_inner_id, x_inner_opt_value in result[X_INNER_OPT].items():
-                row[(X_INNER_OPT, x_inner_id)] = x_inner_opt_value
-
         self._trace = pd.concat(
             (self._trace, pd.DataFrame([row])),
         )

pypesto/objective/amici/amici.py

@@ -14,6 +14,7 @@
     MODE_FUN,
     MODE_RES,
     RDATAS,
+    SPLINE_KNOTS,
     SUFFIXES_CSV,
     SUFFIXES_HDF5,
     ModeType,
@@ -232,8 +233,9 @@ def __init__(
         # `set_custom_timepoints` method for more information.
         self.custom_timepoints = None
 
-        # Initialize the dictionary for saving of inner parameters.
+        # Initialize the list for saving of inner parameter values.
         self.inner_parameters: list[float] = None
+        self.spline_knots: list[list[list[float]]] = None
 
     def get_config(self) -> dict:
         """Return basic information of the objective configuration."""
@@ -504,6 +506,9 @@ def call_unprocessed(
         if ret.get(INNER_PARAMETERS, None) is not None:
             self.inner_parameters = ret[INNER_PARAMETERS]
 
+        if ret.get(SPLINE_KNOTS, None) is not None:
+            self.spline_knots = ret[SPLINE_KNOTS]
+
         # check whether we should update data for preequilibration guesses
         if (
             self.guess_steadystate

pypesto/optimize/optimizer.py

@@ -11,7 +11,7 @@
 import numpy as np
 import scipy.optimize
 
-from ..C import FVAL, GRAD, INNER_PARAMETERS, MODE_FUN, MODE_RES
+from ..C import FVAL, GRAD, INNER_PARAMETERS, MODE_FUN, MODE_RES, SPLINE_KNOTS
 from ..history import HistoryOptions, NoHistory, OptimizerHistory
 from ..objective import Objective
 from ..problem import Problem
@@ -68,6 +68,12 @@ def wrapped_minimize(
         ):
             result[INNER_PARAMETERS] = problem.objective.inner_parameters
 
+        if (
+            hasattr(problem.objective, SPLINE_KNOTS)
+            and problem.objective.spline_knots is not None
+        ):
+            result[SPLINE_KNOTS] = problem.objective.spline_knots
+
         return result
 
     return wrapped_minimize

pypesto/problem/hierarchical.py

@@ -34,6 +34,11 @@ class HierarchicalProblem(Problem):
         Only relevant if hierarchical is True. Contains the bounds of easily
         interpretable inner parameters only, e.g. noise parameters, scaling
         factors, offsets.
+    semiquant_observable_ids:
+        The ids of semiquantitative observables. Only relevant if hierarchical
+        is True. If not None, the optimization result's `spline_knots` will be
+        a list of lists of spline knots for each semiquantitative observable in
+        the order of these ids.
     """
 
     def __init__(
@@ -70,3 +75,7 @@ def __init__(
 
         self.inner_lb = np.array(inner_lb)
         self.inner_ub = np.array(inner_ub)
+
+        self.semiquant_observable_ids = (
+            self.objective.calculator.semiquant_observable_ids
+        )

pypesto/result/optimize.py

@@ -117,6 +117,7 @@ def __init__(
         self.optimizer = optimizer
         self.free_indices = None
         self.inner_parameters = None
+        self.spline_knots = None
 
     def __getattr__(self, key):
         try: