Merge pull request #15 from sijiaww/cvar-upper-only-clean

Add CVaR upper constraint helper in Optimization. Merged to dev branch

Author: gourav3017

portpy/photon/clinical_criteria.py

@@ -271,6 +271,9 @@ def get_dvh_table(self, my_plan: Plan, constraint_list: list = None, opt_params:
         df = pd.DataFrame()
         count = 0
         for i in range(len(dvh_updated_list)):
+
+            dvh_method =  dvh_updated_list[i]['parameters'].get('dvh_method', None)
+
             if 'dose_volume_V' in dvh_updated_list[i]['type']:
                 limit_key = self.matching_keys(dvh_updated_list[i]['constraints'], 'limit')
                 dose_key = self.matching_keys(dvh_updated_list[i]['parameters'], 'dose_')
@@ -279,6 +282,7 @@ def get_dvh_table(self, my_plan: Plan, constraint_list: list = None, opt_params:
                     df.at[count, 'dose_gy'] = self.dose_to_gy(dose_key, dvh_updated_list[i]['parameters'][dose_key])
                     df.at[count, 'volume_perc'] = dvh_updated_list[i]['constraints'][limit_key]
                     df.at[count, 'dvh_type'] = 'constraint'
+                    df.at[count, 'dvh_method'] = dvh_method
                     df.at[count, 'bound_type'] = dvh_updated_list[i]['constraints'].get('bound_type', 'upper')
                     count = count + 1
                 goal_key = self.matching_keys(dvh_updated_list[i]['constraints'], 'goal')
@@ -287,6 +291,7 @@ def get_dvh_table(self, my_plan: Plan, constraint_list: list = None, opt_params:
                     df.at[count, 'dose_gy'] = self.dose_to_gy(dose_key, dvh_updated_list[i]['parameters'][dose_key])
                     df.at[count, 'volume_perc'] = dvh_updated_list[i]['constraints'][goal_key]
                     df.at[count, 'dvh_type'] = 'goal'
+                    df.at[count, 'dvh_method'] = dvh_method
                     df.at[count, 'weight'] = dvh_updated_list[i]['parameters']['weight']
                     df.at[count, 'bound_type'] = dvh_updated_list[i]['constraints'].get('bound_type', 'upper')
                     count = count + 1
@@ -296,6 +301,7 @@ def get_dvh_table(self, my_plan: Plan, constraint_list: list = None, opt_params:
                     df.at[count, 'structure_name'] = dvh_updated_list[i]['parameters']['structure_name']
                     df.at[count, 'volume_perc'] = dvh_updated_list[i]['parameters']['volume_perc']
                     df.at[count, 'dose_gy'] = self.dose_to_gy(limit_key, dvh_updated_list[i]['constraints'][limit_key])
+                    df.at[count, 'dvh_method'] = dvh_method
                     df.at[count, 'dvh_type'] = 'constraint'
                     df.at[count, 'bound_type'] = dvh_updated_list[i]['constraints'].get('bound_type', 'upper')
                     count = count + 1
@@ -304,6 +310,7 @@ def get_dvh_table(self, my_plan: Plan, constraint_list: list = None, opt_params:
                     df.at[count, 'structure_name'] = dvh_updated_list[i]['parameters']['structure_name']
                     df.at[count, 'volume_perc'] = dvh_updated_list[i]['parameters']['volume_perc']
                     df.at[count, 'dose_gy'] = self.dose_to_gy(goal_key, dvh_updated_list[i]['constraints'][goal_key])
+                    df.at[count, 'dvh_method'] = dvh_method
                     df.at[count, 'dvh_type'] = 'goal'
                     df.at[count, 'weight'] = dvh_updated_list[i]['parameters']['weight']
                     df.at[count, 'bound_type'] = dvh_updated_list[i]['constraints'].get('bound_type', 'upper')

portpy/photon/optimization.py

@@ -159,8 +159,73 @@ def create_cvxpy_problem(self):
         d_max = np.infty * np.ones(A.shape[0]) # create a vector to avoid putting redundant max constraint on
         # duplicate voxels among structure
 
-        # Adding max/mean constraints
+        # Adding constraints
         for i in range(len(constraint_def)):
+
+            item = constraint_def[i]
+            item_type = item.get('type')
+
+            if item_type in ('dose_volume_V', 'dose_volume_D'):
+
+                params = item.get('parameters', {})
+                cons   = item.get('constraints', {})
+
+                dvh_method = params.get('dvh_method', None)
+                if dvh_method != 'cvar':
+                    continue
+
+                struct = params['structure_name']
+                voxel_idx = st.get_opt_voxels_idx(struct)
+                if len(voxel_idx) == 0:
+                    continue
+
+                limit = float(params['dose_gy'])
+
+                volume_perc = float(cons['limit_volume_perc'])
+                if not (0.0 < volume_perc < 100.0):
+                    raise ValueError(
+                        f"limit_volume_perc must be in (0,100), got {volume_perc}"
+                    )
+
+                alpha = 1.0 - volume_perc / 100.0
+                dose_1d_list = A[voxel_idx, :] @ x * num_fractions
+
+                if cons.get('constraint_type') == 'upper':
+
+                    label = f"constraint_{struct}_{item_type}_upper_a{alpha:.4f}"
+
+                    zeta = cp.Variable(name=f"zeta_{label}")
+                    w = cp.Variable(len(voxel_idx), name=f"w_{label}")
+
+                    self.vars[f"zeta_{label}"] = zeta
+                    self.vars[f"w_{label}"] = w
+
+                    constraints += [
+                        zeta + (1.0 / ((1.0 - alpha) * len(voxel_idx))) * cp.sum(w) <= limit,
+                        w >= dose_1d_list - zeta,
+                        w >= 0
+                    ]
+
+                elif cons.get('constraint_type') == 'lower':
+
+                    dose_neg = -dose_1d_list
+                    limit_neg = -limit
+
+                    label = f"constraint_{struct}_{item_type}_lower_a{alpha:.4f}"
+
+                    zeta = cp.Variable(name=f"zeta_{label}")
+                    w = cp.Variable(len(voxel_idx), name=f"w_{label}")
+
+                    self.vars[f"zeta_{label}"] = zeta
+                    self.vars[f"w_{label}"] = w
+
+                    constraints += [
+                        zeta + (1.0 / ((1.0 - alpha) * len(voxel_idx))) * cp.sum(w) <= limit_neg,
+                        w >= dose_neg - zeta,
+                        w >= 0
+                    ]
+
+
             if constraint_def[i]['type'] == 'max_dose':
                 limit_key = self.matching_keys(constraint_def[i]['constraints'], 'limit')
                 if limit_key:
@@ -188,13 +253,15 @@ def create_cvxpy_problem(self):
                                         (cp.sum((cp.multiply(st.get_opt_voxels_volume_cc(org),
                                                              A[st.get_opt_voxels_idx(org), :] @ x))))
                                         <= limit / num_fractions]
+
         mask = np.isfinite(d_max)
         # Create index mask arrays
         indices = np.arange(len(mask))  # assumes mask is 1D and corresponds to voxel indices
         all_d_max_vox_ind = indices[mask]
         constraints += [A[all_d_max_vox_ind, :] @ x <= d_max[all_d_max_vox_ind]] # Add constraint for all d_max voxels at once
         print('Problem created')
 
+
     def add_max(self, struct: str, dose_gy: float):
         """
         Add max constraints to the problem