More flexible optimization framework

pyRadPlan/optimization/new_templates/objectives.py

@@ -0,0 +1,21 @@
+import numpy as np
+
+
+class abstract_objective:
+    def __init__(self):
+        pass
+
+    def evaluate(x: np.ndarray[float]) -> float:
+        pass
+
+    def evaluate_gradient(x: np.ndarray[float]) -> np.ndarray[float]:
+        pass
+
+    def evaluate_hessian(x: np.ndarray[float]) -> np.ndarray[float]:
+        pass
+
+    def is_linear() -> bool:
+        pass
+
+    def is_convex() -> bool:
+        pass

pyRadPlan/optimization/new_templates/planning_problem_sketch.py

@@ -0,0 +1,53 @@
+import numpy as np
+from enum import Enum
+
+
+class PlanningProblem:
+    def __init__(
+        self,
+        scalarization_method: str | Enum,
+        tradeoff_exploration_method: str | Enum | None,
+        method_parameters: dict,
+    ):
+        pass
+
+    def evaluate_objectives(x: np.ndarray[float]) -> list[np.ndarray[float]]:
+        pass
+
+    def evaluate_objective_gradients(x: np.ndarray[float]) -> list[np.ndarray[float]]:
+        # returns list of 2d arrays
+        pass
+
+    def evaluate_objective_hessian(x: np.ndarray[float]) -> list[np.ndarray[float]]:
+        # returns list of 3d arrays
+        pass
+
+    def evaluate_constraints(x: np.ndarray[float]) -> list[np.ndarray[float]]:
+        pass
+
+    def evaluate_constraints_gradients(x: np.ndarray[float]) -> list[np.ndarray[float]]:
+        # returns list of 2d arrays
+        pass
+
+    def evaluate_constraints_hessian(x: np.ndarray[float]) -> list[np.ndarray[float]]:
+        # returns list of 3d arrays
+        pass
+
+    def are_objectives_convex() -> bool:
+        pass
+
+    def are_objectives_linear() -> bool:
+        pass
+
+    def solve(y: np.ndarray[float]) -> list[np.ndarray[float]]:
+        pass
+
+    def make_tradeoff_exploration_method_instance(
+        evaluate_objectives: callable[np.ndarray[float], list[np.ndarray[float]]],
+        evaluate_constraints: callable[np.ndarray[float], list[np.ndarray[float]]],
+        evaluate_x_gradients,
+        etc,
+        scalarization_method: str,
+        method_parameters: dict,
+    ) -> TradeoffExplorationMethod:
+        pass

pyRadPlan/optimization/new_templates/scalarization_method.py

@@ -0,0 +1,44 @@
+import numpy as np
+
+
+class ScalarizationMethod:
+    def __init__(
+        self,
+        evaluate_objectives: callable[np.ndarray[float], list[np.ndarray[float]]],
+        evaluate_constraints: callable[np.ndarray[float], list[np.ndarray[float]]],
+        evaluate_x_gradients,
+        etc,
+        parameters: dict,
+    ):
+        # Implementations of class should also manage the concatenating the additional variables and separating them
+        pass
+
+    def variable_upper_bounds() -> np.ndarray[float]:
+        pass
+
+    def variable_lower_bounds() -> np.ndarray[float]:
+        pass
+
+    def get_linear_constrains(self) -> dict[Index, LinearConstraint]:
+        pass
+
+    def get_nonlinear_constraints(self) -> dict[Index, NonlinearConstraint]:
+        pass
+
+    def evaluate_objective(x: np.ndarray) -> float:
+        print("This is not the same as evaluating objectives. E.g. make weighted sum")
+
+    def evaluate_constraints(x: np.ndarray) -> np.ndarray:
+        print(
+            "Most of the time this will be the objective constraints and the constraints from the scalarization method"
+        )
+
+    def solve(self, x: np.ndarray[float]) -> np.ndarray[float]:
+        print("Do something")
+        self._call_solver_interface(self.solver, self.solver_params)
+
+    def is_objective_convex() -> bool:
+        pass
+
+    def _call_solver_interface(solver: str, params: dict) -> np.ndarray[float]:
+        pass

pyRadPlan/optimization/new_templates/tradeoff_exploration_method.py

@@ -0,0 +1,9 @@
+import numpy as np
+
+
+class TradeoffExplorationMethod:
+    def __init__(self, params: dict):
+        pass
+
+    def solve(x: np.ndarray[float]) -> list[np.ndarray[float]]:
+        pass

pyRadPlan/optimization/objectives/_max_dvh.py

@@ -22,7 +22,8 @@ class MaxDVH(Objective):
     """
 
     name = "Max DVH"
-
+    is_convex = False
+    is_linear = True
     d: Annotated[float, Field(default=30.0, ge=0.0), ParameterMetadata(kind="reference")]
     v_max: Annotated[
         float, Field(default=50.0, ge=0.0, le=100.0), ParameterMetadata(kind="relative_volume")

pyRadPlan/optimization/objectives/_min_dvh.py

@@ -22,7 +22,8 @@ class MinDVH(Objective):
     """
 
     name = "Min DVH"
-
+    is_convex = False
+    is_linear = True
     d: Annotated[float, Field(default=30.0, ge=0.0), ParameterMetadata(kind="reference")]
     v_min: Annotated[
         float, Field(default=95.0, ge=0.0, le=100.0), ParameterMetadata(kind="relative_volume")

pyRadPlan/optimization/objectives/_objective.py

@@ -72,10 +72,18 @@ class Objective(PyRadPlanBaseModel):
         Weight/Priority assigned to the objective function.
     quantity : str
         The quantity this objective is connected to (e.g. 'physical_dose', 'RBExDose').
+    is_linear : bool
+        bool to indicate if the objective is linear.
+    is_convex : bool
+        bool to indicate if the objective is convex.
     """
 
     name: ClassVar[str]
     has_hessian: ClassVar[bool] = False
+
+    is_linear: ClassVar[bool] = False
+    is_convex: ClassVar[bool] = True
+
     priority: float = Field(default=1.0, ge=0.0, alias="penalty")
     quantity: str = Field(default="physical_dose")
 

pyRadPlan/optimization/problems/_nonlin_fluence.py

@@ -8,7 +8,6 @@
 
 from ...plan import Plan
 from ._optiprob import NonLinearPlanningProblem
-from ..solvers import NonLinearOptimizer
 from ..objectives import Objective
 
 logger = logging.getLogger(__name__)
@@ -37,15 +36,13 @@ class NonLinearFluencePlanningProblem(NonLinearPlanningProblem):
     bypass_objective_jacobian: bool
 
     def __init__(self, pln: Union[Plan, dict] = None):
-        self.bypass_objective_jacobian = True
+        self.bypass_objective_jacobian = False
 
         self._target_voxels = None
         self._patient_voxels = None
 
         self._grad_cache_intermediate = None
         self._grad_cache = None
-        self._obj_times = []
-        self._deriv_times = []
         self._solve_time = None
 
         super().__init__(pln)
@@ -54,17 +51,17 @@ def _initialize(self):
         """Initialize this problem."""
         super()._initialize()
 
-        # Check if the solver is adequate to solve this problem
-        # TODO: check that it can do constraints
-        if not isinstance(self.solver, NonLinearOptimizer):
-            raise ValueError("Solver must be an instance of SolverBase")
-
-        self.solver.objective = self._objective_function
-        self.solver.gradient = self._objective_gradient
-        self.solver.bounds = (0.0, np.inf)
-        self.solver.max_iter = 500
-
-    def _objective_functions(self, x: NDArray) -> NDArray:
+        # # Check if the solver is adequate to solve this problem
+        # # TODO: check that it can do constraints
+        # # if not isinstance(self.solver, NonLinearOptimizer):
+        # #     raise ValueError("Solver must be an instance of SolverBase")
+        # self.solver = get_solver(self.solver)
+        # self.solver.objective = self._objective_function
+        # self.solver.gradient = self._objective_gradient
+        # self.solver.bounds = (0.0, np.inf)
+        # self.solver.max_iter = 500
+
+    def _evaluate_objective_functions(self, x: NDArray) -> NDArray:
         """Define the objective functions."""
 
         q_vectors = {}
@@ -85,8 +82,7 @@ def _objective_functions(self, x: NDArray) -> NDArray:
                 f_vals.append(
                     sum(
                         [
-                            obj.priority
-                            * obj.compute_objective(q_vectors[obj.quantity].flat[scen_ix][ix])
+                            obj.compute_objective(q_vectors[obj.quantity].flat[scen_ix][ix])
                             for scen_ix in q_scenarios[obj.quantity]
                         ]
                     )
@@ -95,13 +91,11 @@ def _objective_functions(self, x: NDArray) -> NDArray:
         # return as numpy array
         return np.asarray(f_vals, dtype=np.float64)
 
-    def _objective_function(self, x: NDArray) -> np.float64:
-        t = time.time()
-        f = np.sum(self._objective_functions(x))
-        self._obj_times.append(time.time() - t)
-        return f
+    # def _objective_function(self, x: NDArray) -> np.float64:
+    #     f = np.sum(self._evaluate_objective_functions(x))
+    #     return f
 
-    def _objective_jacobian(self, x: NDArray) -> NDArray:
+    def _evaluate_objective_jacobian(self, x: NDArray) -> NDArray:
         """Define the objective jacobian."""
 
         q_vectors = {}
@@ -142,8 +136,7 @@ def _objective_jacobian(self, x: NDArray) -> NDArray:
                     q_cache_index = self._q_cache_index[cnt]
                 for scen_ix in q_scenarios[obj.quantity]:
                     self._grad_cache_intermediate[obj.quantity][q_cache_index, ix] += (
-                        obj.priority
-                        * obj.compute_gradient(q_vectors[obj.quantity].flat[scen_ix][ix])
+                        obj.compute_gradient(q_vectors[obj.quantity].flat[scen_ix][ix])
                     )
                 cnt += 1
 
@@ -175,55 +168,35 @@ def _objective_jacobian(self, x: NDArray) -> NDArray:
 
         return self._grad_cache
 
-    def _objective_gradient(self, x: NDArray) -> NDArray:
-        t = time.time()
-        jac = np.sum(self._objective_jacobian(x), axis=0)
-        self._deriv_times.append(time.time() - t)
-        return jac
-
-    def _objective_hessian(self, x: NDArray) -> NDArray:
+    def _evaluate_objective_hessian(self, x: NDArray) -> NDArray:
         """Define the objective hessian."""
         return {}
 
-    def _constraint_functions(self, x: NDArray) -> NDArray:
+    def _evaluate_constraint_functions(self, x: NDArray) -> NDArray:
         """Define the constraint functions."""
         return None
 
-    def _constraint_jacobian(self, x: NDArray) -> NDArray:
+    def _evaluate_constraint_jacobian(self, x: NDArray) -> NDArray:
         """Define the constraint jacobian."""
         return None
 
-    def _constraint_jacobian_structure(self) -> NDArray:
+    def _evaluate_constraint_jacobian_structure(self) -> NDArray:
         """Define the constraint jacobian structure."""
         return None
 
-    def _variable_bounds(self, x: NDArray) -> NDArray:
+    def _get_variable_bounds(self, x: NDArray) -> NDArray:
         """Define the variable bounds."""
         return {}
 
     def _solve(self) -> tuple[NDArray, dict]:
         """Solve the problem."""
 
-        self._deriv_times = []
-        self._obj_times = []
-
         x0 = np.zeros((self._dij.total_num_of_bixels,), dtype=np.float64)
         t = time.time()
-        result = self.solver.solve(x0)
+        result = self.tradeoff_strategy.solve(x0)
+        # result = self.solver.solve(x0)
         self._solve_time = time.time() - t
 
-        logger.info(
-            "%d Objective function evaluations, avg. time: %g +/- %g s",
-            len(self._obj_times),
-            np.mean(self._obj_times),
-            np.std(self._obj_times),
-        )
-        logger.info(
-            "%d Derivative evaluations, avg. time: %g +/- %g s",
-            len(self._deriv_times),
-            np.mean(self._deriv_times),
-            np.std(self._deriv_times),
-        )
         logger.info("Solver time: %g s", self._solve_time)
 
         return result