diff --git a/portpy/photon/clinical_criteria.py b/portpy/photon/clinical_criteria.py index 6f1d2d5..27940c5 100644 --- a/portpy/photon/clinical_criteria.py +++ b/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') diff --git a/portpy/photon/optimization.py b/portpy/photon/optimization.py index aaeff72..8e81dd6 100644 --- a/portpy/photon/optimization.py +++ b/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,6 +253,7 @@ 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 @@ -195,6 +261,7 @@ def create_cvxpy_problem(self): 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