Fix 1474 - add fields to result object - correct calculation of number of clps.

glotaran/optimization/objective.py

@@ -1,5 +1,6 @@
 from __future__ import annotations
 
+from dataclasses import dataclass
 from typing import TYPE_CHECKING
 
 import numpy as np
@@ -19,6 +20,14 @@
     from glotaran.typing.types import ArrayLike
 
 
+@dataclass
+class OptimizationObjectiveResult:
+    data: dict[str, xr.Dataset]
+    free_clp_size: int
+    additional_penalty: float
+    dataset_penalty: dict[str, float]
+
+
 class OptimizationObjective:
     def __init__(self, model: ExperimentModel):
         self._data = (
@@ -90,7 +99,7 @@ def calculate(self) -> ArrayLike:
             )
         return np.concatenate(penalties)
 
-    def get_result(self) -> dict[str, xr.Dataset]:
+    def get_result(self) -> OptimizationObjectiveResult:
         return (
             self.create_unlinked_result()
             if isinstance(self._data, OptimizationData)
@@ -268,14 +277,15 @@ def finalize_result_dataset(self, dataset: xr.Dataset, data: OptimizationData, a
                 data.model, dataset, True
             )
 
-    def create_linked_result(self) -> dict[str, xr.Dataset]:
+    def create_linked_result(self) -> OptimizationObjectiveResult:
         assert isinstance(self._data, LinkedOptimizationData)
         matrices = {
             label: OptimizationMatrix.from_data(data) for label, data in self._data.data.items()
         }
         linked_matrices = OptimizationMatrix.from_linked_data(self._data, matrices)
         clp_axes = [matrix.clp_axis for matrix in linked_matrices]
         reduced_matrices = self.calculate_reduced_matrices(linked_matrices)
+        free_clp_size = sum(len(matrix.clp_axis) for matrix in reduced_matrices)
         estimations = self.resolve_estimations(
             linked_matrices,
             reduced_matrices,
@@ -291,28 +301,54 @@ def create_linked_result(self) -> dict[str, xr.Dataset]:
                 results[label], label, clp_axes, estimations
             )
             self.finalize_result_dataset(results[label], data)
-        return results
+        additional_penalty = sum(
+            calculate_clp_penalties(
+                linked_matrices,
+                estimations,
+                self._data.global_axis,
+                self._model.clp_penalties,
+            )
+        )
+        dataset_penalties = {label: d.residual.sum().data for label, d in results.items()}
+        return OptimizationObjectiveResult(
+            results, free_clp_size, additional_penalty, dataset_penalties
+        )
 
-    def create_unlinked_result(self) -> dict[str, xr.Dataset]:
+    def create_unlinked_result(self) -> OptimizationObjectiveResult:
         assert isinstance(self._data, OptimizationData)
 
         label = next(iter(self._model.datasets.keys()))
         result = self.create_result_dataset(label, self._data)
 
         matrix = OptimizationMatrix.from_data(self._data)
         self.add_matrix_to_dataset(result, matrix)
+        additional_penalty = 0
         if self._data.is_global:
             self.add_global_clp_and_residual_to_dataset(result, self._data, matrix)
+            free_clp_size = len(matrix.clp_axis)
+
         else:
             reduced_matrices = self.calculate_reduced_matrices(
                 matrix.as_global_list(self._data.global_axis)
             )
+            free_clp_size = sum(len(matrix.clp_axis) for matrix in reduced_matrices)
             estimations = self.resolve_estimations(
                 matrix.as_global_list(self._data.global_axis),
                 reduced_matrices,
                 self.calculate_estimations(reduced_matrices),
             )
             self.add_unlinked_clp_and_residual_to_dataset(result, estimations)
+            additional_penalty = sum(
+                calculate_clp_penalties(
+                    [matrix],
+                    estimations,
+                    self._data.global_axis,
+                    self._model.clp_penalties,
+                )
+            )
 
         self.finalize_result_dataset(result, self._data)
-        return {label: result}
+        dataset_penalties = {label: result.residual.sum().data}
+        return OptimizationObjectiveResult(
+            {label: result}, free_clp_size, additional_penalty, dataset_penalties
+        )

glotaran/optimization/optimization.py

@@ -135,33 +135,39 @@ def run(self) -> tuple[Parameters, dict[str, xr.Dataset], OptimizationResult]:
                 termination_reason = str(e)
 
         penalty = np.concatenate([o.calculate() for o in self._objectives])
-        data = dict(ChainMap(*[o.get_result() for o in self._objectives]))
-        nr_clp = len({str(c.data) for d in data.values() for c in d.clp_label})
+        results = [o.get_result() for o in self._objectives]
+        data = dict(ChainMap(*[r.data for r in results]))
+        number_of_clps = sum(r.free_clp_size for r in results)
+        additional_penalty = sum(r.additional_penalty for r in results)
         result = OptimizationResult.from_least_squares_result(
             ls_result,
             self._parameter_history,
             OptimizationHistory.from_stdout_str(self._tee.read()),
             penalty,
+            additional_penalty,
             self._free_parameter_labels,
             termination_reason,
-            nr_clp,
+            number_of_clps,
         )
         return self._parameters, data, result
 
     def dry_run(self) -> tuple[Parameters, dict[str, xr.Dataset], OptimizationResult]:
         termination_reason = "Dry run."
 
         penalty = np.concatenate([o.calculate() for o in self._objectives])
-        data = dict(ChainMap(*[o.get_result() for o in self._objectives]))
-        nr_clp = len({str(c.data) for d in data.values() for c in d.clp_label})
+        results = [o.get_result() for o in self._objectives]
+        data = dict(ChainMap(*[r.data for r in results]))
+        number_of_clps = sum(r.free_clp_size for r in results)
+        additional_penalty = sum(r.additional_penalty for r in results)
         result = OptimizationResult.from_least_squares_result(
             None,
             self._parameter_history,
             OptimizationHistory.from_stdout_str(self._tee.read()),
             penalty,
+            additional_penalty,
             self._free_parameter_labels,
             termination_reason,
-            nr_clp,
+            number_of_clps,
         )
         return self._parameters, data, result
 

glotaran/optimization/result.py

@@ -61,7 +61,7 @@ class OptimizationResult(BaseModel):
 
     The rows and columns are corresponding to :attr:`free_parameter_labels`."""
 
-    number_clp: int | None = None
+    number_of_clps: int | None = None
     degrees_of_freedom: int | None = None
     """Degrees of freedom in optimization :math:`N - N_{vars}`."""
 
@@ -88,6 +88,7 @@ class OptimizationResult(BaseModel):
 
     :math:`rms = \sqrt{\chi^2_{red}}`
     """
+    additional_penalty: float | None = None
 
     @classmethod
     def from_least_squares_result(
@@ -96,9 +97,10 @@ def from_least_squares_result(
         parameter_history: ParameterHistory,
         optimization_history: OptimizationHistory,
         penalty: ArrayLike,
+        additional_penalty: float,
         free_parameter_labels: list[str],
         termination_reason: str,
-        number_clp: int,
+        number_of_clps: int,
     ):
         success = result is not None
 
@@ -115,15 +117,16 @@ def from_least_squares_result(
         }
 
         if success:
-            result_args["number_clp"] = number_clp
+            result_args["number_of_clps"] = number_of_clps
+            result_args["additional_penalty"] = additional_penalty
             result_args["number_of_jacobian_evaluations"] = result.njev  # type:ignore[union-attr]
             result_args["optimality"] = float(result.optimality)  # type:ignore[union-attr]
             result_args["number_of_data_points"] = result.fun.size  # type:ignore[union-attr]
             result_args["number_of_parameters"] = result.x.size  # type:ignore[union-attr]
             result_args["degrees_of_freedom"] = (
                 result_args["number_of_data_points"]
                 - result_args["number_of_parameters"]
-                - result_args["number_clp"]
+                - result_args["number_of_clps"]
             )
             result_args["chi_square"] = float(np.sum(result.fun**2))  # type:ignore[union-attr]
             result_args["reduced_chi_square"] = (

tests/optimization/test_objective.py

@@ -26,7 +26,7 @@ def test_single_data():
     data_size = data_model.data["model"].size * data_model.data["global"].size
     assert penalty.size == data_size
 
-    result = objective.get_result()
+    result = objective.get_result().data
     assert "test" in result
 
     result_data = result["test"]
@@ -51,7 +51,7 @@ def test_global_data(weight: bool):
     data_size = data_model.data["model"].size * data_model.data["global"].size
     assert penalty.size == data_size
 
-    result = objective.get_result()
+    result = objective.get_result().data
     assert "test" in result
 
     result_data = result["test"]
@@ -90,7 +90,7 @@ def test_multiple_data():
     data_size_two = data_model_two.data["model"].size * data_model_two.data["global"].size
     assert penalty.size == data_size_one + data_size_two
 
-    result = objective.get_result()
+    result = objective.get_result().data
 
     assert "independent" in result
     result_data = result["independent"]
@@ -129,7 +129,7 @@ def test_result_data(weight: bool):
     data_size = data_model.data["model"].size * data_model.data["global"].size
     assert penalty.size == data_size
 
-    result = objective.get_result()
+    result = objective.get_result().data
     assert "test" in result
 
     result_data = result["test"]