Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Update to store objectives and constraints as optimizationProblem properties #648

Closed
wants to merge 5 commits into from
Closed

Update to store objectives and constraints as optimizationProblem properties #648

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
fe06707
Updates to where objectives and constraints are stored
tobiasbecher May 31, 2023
823794f
Merge branch 'dev_Pareto' into dev_ParetoObj
tobiasbecher Jun 5, 2023
4cba6b2
Merge branch 'dev_Pareto' into dev_ParetoObj
tobiasbecher Jul 21, 2023
2e806a9
Merge branch 'dev_varRBErobOpt' into dev_ParetoObj
tobiasbecher Jul 21, 2023
bc7f21d
Fixed bug with constraints and added some comments
tobiasbecher Jul 28, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 2 additions & 1 deletion matRad_fluenceOptimization.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -249,7 +249,7 @@
backProjection.scenarios = ixForOpt;
backProjection.scenarioProb = pln.multScen.scenProb;

optiProb = matRad_OptimizationProblem(backProjection);
optiProb = matRad_OptimizationProblem(backProjection,cst);
tobiasbecher marked this conversation as resolved.
Show resolved Hide resolved
optiProb.quantityOpt = pln.bioParam.quantityOpt;
if isfield(pln,'propOpt') && isfield(pln.propOpt,'useLogSumExpForRobOpt')
optiProb.useLogSumExpForRobOpt = pln.propOpt.useLogSumExpForRobOpt;
Expand Down Expand Up @@ -307,6 +307,7 @@
resultGUI.wUnsequenced = wOpt;
resultGUI.usedOptimizer = optimizer;
resultGUI.info = info;
resultGUI.optiProb = optiProb;

%Robust quantities
if FLAG_ROB_OPT || numel(ixForOpt) > 1
Expand Down
98 changes: 94 additions & 4 deletions optimization/@matRad_OptimizationProblem/matRad_OptimizationProblem.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,17 +24,27 @@

properties
BP
normalizationScheme = struct('type','none');
tobiasbecher marked this conversation as resolved.
Show resolved Hide resolved
objectives = {}; %cell array storing all objectives, has to be initialized at the start
constraints = {}; %
objidx;
tobiasbecher marked this conversation as resolved.
Show resolved Hide resolved
constridx;
quantityOpt = '';
useMaxApprox = 'logsumexp'; %'pnorm'; %'logsumexp'; %'none';
p = 30; %Can be chosen larger (closer to maximum) or smaller (closer to mean). Only tested 20 >= p >= 1

minimumW = NaN;
maximumW = NaN;
end


methods
function obj = matRad_OptimizationProblem(backProjection)
obj.BP = backProjection;
function obj = matRad_OptimizationProblem(backProjection,cst)
tobiasbecher marked this conversation as resolved.
Show resolved Hide resolved

obj.BP = backProjection; %needs to be initalized to have access to setBiologicalDosePrescriptions
if nargin == 2
obj.extractObjectivesAndConstraintsFromcst(cst);
end
end

%Objective function declaration
Expand Down Expand Up @@ -81,6 +91,50 @@
matRad_cfg.dispError('Maximum Bounds for Optimization Problem could not be set!');
end
end

function normalizedfVals = normalizeObjectives(optiProb,fVals)
%function to normalize objectives (used for sandwich
%algorithms)
switch optiProb.normalizationScheme.type
case 'none'
%default case no normalization
normalizedfVals = fVals;
case 'UL'
%used to normalize with respect to min and max values
%maybe check that U and L are defined
normalizedfVals = (fVals - optiProb.normalizationScheme.L)./(optiProb.normalizationScheme.U-optiProb.normalizationScheme.L); %might have to check that U and L work!
otherwise
matRad_cfg.dispError('Normalization scheme not known!');
end
end

function normalizedGradient = normalizeGradients(optiProb,Gradient)
%function to normalize objectives (used for sandwich
%algorithms)
switch optiProb.normalizationScheme.type

case 'none'
%default case no normalization
normalizedGradient = Gradient;
case 'UL'
%used to normalize with respect to min and max values
%maybe check that U and L are defined
normalizedGradient = Gradient./(optiProb.normalizationScheme.U-optiProb.normalizationScheme.L); %might have to check that U and L work!
otherwise
matRad_cfg.dispError('Normalization scheme not known!');
end
end

function updatePenalties(optiProb,newPen)
%TODO: Allow grouping of objectives
if numel(optiProb.objectives) ~= numel(newPen)
matRad_cfg.dispError('Number of objectives in optimization Problem not equal to number of new penalties to be set!');
end
for i=1:numel(newPen)
optiProb.objectives{i}.penalty = newPen(i)*100; %scaling of penalties with factor 100
end
end

end

methods (Access = protected)
Expand Down Expand Up @@ -117,6 +171,42 @@

grad = fac * (tmp ./ pNormVal).^(p-1);
end
end
end
end


methods (Access = private)
function extractObjectivesAndConstraintsFromcst(optiProb,cst)
%used to extract objectives from cst and store in cell array as property of optimization Problem
optiProb.objidx = [];
optiProb.constridx = [];
optiProb.objectives = {};
optiProb.constraints = {};

for i = 1:size(cst,1) % loop over cst

if ~isempty(cst{i,4}{1}) && ( isequal(cst{i,3},'OAR') || isequal(cst{i,3},'TARGET') )

for j = 1:numel(cst{i,6})
%check whether dose objective or constraint
obj = cst{i,6}{j};
if isstruct(cst{i,6}{j})
obj = matRad_DoseOptimizationFunction.createInstanceFromStruct(obj);
end
if contains(class(obj),'DoseObjectives')
optiProb.objidx = [optiProb.objidx;i,j];
obj = optiProb.BP.setBiologicalDosePrescriptions(obj,cst{i,5}.alphaX,cst{i,5}.betaX);
optiProb.objectives(end+1) = {obj};

elseif contains(class(obj),'DoseConstraints')
optiProb.constridx = [optiProb.constridx;i,j];
obj = optiProb.BP.setBiologicalDosePrescriptions(obj,cst{i,5}.alphaX,cst{i,5}.betaX);
optiProb.constraints(end+1) = {obj};
end
end
end
end
end


end
end
150 changes: 65 additions & 85 deletions optimization/@matRad_OptimizationProblem/matRad_constraintFunctions.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,6 @@
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% get current dose / effect / RBExDose vector
optiProb.BP.compute(dij,w);
d = optiProb.BP.GetResult();
Expand All @@ -50,96 +49,77 @@
c = [];

% compute objective function for every VOI.
for i = 1:size(cst,1)

% Only take OAR or target VOI.
if ~isempty(cst{i,4}{1}) && ( isequal(cst{i,3},'OAR') || isequal(cst{i,3},'TARGET') )

% loop over the number of constraints for the current VOI
for j = 1:numel(cst{i,6})
for i = 1:size(optiProb.constridx,1)

constraint = cst{i,6}{j};
constraint = optiProb.constraints{i};
curConidx = optiProb.constridx(i,1);

robustness = constraint.robustness;

% only perform computations for constraints
if isa(constraint,'DoseConstraints.matRad_DoseConstraint')
switch robustness
case 'none' % if conventional opt: just sum objectives of nominal dose
d_i = d{1}(cst{curConidx,4}{1});
c = [c; constraint.computeDoseConstraintFunction(d_i)];

case 'PROB' % if prob opt: sum up expectation value of objectives

d_i = dExp{1}(cst{curConidx,4}{1});
c = [c; constraint.computeDoseConstraintFunction(d_i)];

case 'VWWC' % voxel-wise worst case - takes minimum dose in TARGET and maximum in OAR
contourIx = unique(contourScen);
if ~isscalar(contourIx)
% voxels need to be tracked through the 4D CT,
% not yet implemented
matRad_cfg.dispError('4D VWWC optimization is currently not supported');
end

% prepare min/max dose vector
if ~exist('d_tmp','var')
d_tmp = [d{useScen}];
end

% rescale dose parameters to biological optimization quantity if required
constraint = optiProb.BP.setBiologicalDosePrescriptions(constraint,cst{i,5}.alphaX,cst{i,5}.betaX);
d_Scen = d_tmp(cst{curConidx,4}{contourIx},:);

% retrieve the robustness type
robustness = constraint.robustness;
d_max = max(d_Scen,[],2);
d_min = min(d_Scen,[],2);

switch robustness
case 'none' % if conventional opt: just sum objectives of nominal dose
d_i = d{1}(cst{i,4}{1});
c = [c; constraint.computeDoseConstraintFunction(d_i)];

case 'PROB' % if prob opt: sum up expectation value of objectives

d_i = dExp{1}(cst{i,4}{1});
c = [c; constraint.computeDoseConstraintFunction(d_i)];

case 'VWWC' % voxel-wise worst case - takes minimum dose in TARGET and maximum in OAR
contourIx = unique(contourScen);
if ~isscalar(contourIx)
% voxels need to be tracked through the 4D CT,
% not yet implemented
matRad_cfg.dispError('4D VWWC optimization is currently not supported');
end

% prepare min/max dose vector
if ~exist('d_tmp','var')
d_tmp = [d{useScen}];
end

d_Scen = d_tmp(cst{i,4}{contourIx},:);

d_max = max(d_Scen,[],2);
d_min = min(d_Scen,[],2);

if isequal(cst{i,3},'OAR')
d_i = d_max;
elseif isequal(cst{i,3},'TARGET')
d_i = d_min;
end

c = [c; constraint.computeDoseConstraintFunction(d_i)];

case 'VWWC_INV' %inverse voxel-wise conformitiy - takes maximum dose in TARGET and minimum in OAR
contourIx = unique(contourScen);
if ~isscalar(contourIx)
% voxels need to be tracked through the 4D CT,
% not yet implemented
matRad_cfg.dispError('4D inverted VWWC optimization is currently not supported');
end

% prepare min/max dose vector
if ~exist('d_tmp','var')
d_tmp = [d{:}];
end

d_Scen = d_tmp(cst{i,4}{contourIx},:);
d_max = max(d_Scen,[],2);
d_min = min(d_Scen,[],2);

if isequal(cst{i,3},'OAR')
d_i = d_min;
elseif isequal(cst{i,3},'TARGET')
d_i = d_max;
end

c = [c; constraint.computeDoseConstraintFunction(d_i)];
otherwise
matRad_cfg.dispError('Robustness setting %s not yet supported!',constraint.robustness);
if isequal(cst{curConidx,3},'OAR')
d_i = d_max;
elseif isequal(cst{curConidx,3},'TARGET')
d_i = d_min;
end

c = [c; constraint.computeDoseConstraintFunction(d_i)];

end
case 'VWWC_INV' %inverse voxel-wise conformitiy - takes maximum dose in TARGET and minimum in OAR
contourIx = unique(contourScen);
if ~isscalar(contourIx)
% voxels need to be tracked through the 4D CT,
% not yet implemented
matRad_cfg.dispError('4D inverted VWWC optimization is currently not supported');
end

% prepare min/max dose vector
if ~exist('d_tmp','var')
d_tmp = [d{:}];
end

d_Scen = d_tmp(cst{curConidx,4}{contourIx},:);
d_max = max(d_Scen,[],2);
d_min = min(d_Scen,[],2);

if isequal(cst{curConidx,3},'OAR')
d_i = d_min;
elseif isequal(cst{curConidx,3},'TARGET')
d_i = d_max;
end

c = [c; constraint.computeDoseConstraintFunction(d_i)];
otherwise
matRad_cfg.dispError('Robustness setting %s not yet supported!',constraint.robustness);
end

end % if we are a constraint

end % over all defined constraints & objectives

end % if structure not empty and oar or target

end % over all structures
end % if we are a constraint
end
% if structure not empty and oar or target
Loading
Loading