matRad_MinMaxDose

0001 classdef matRad_MinMaxDose < DoseConstraints.matRad_DoseConstraint
0002     % matRad_MinMaxDose Implements a MinMaxDose constraint
0003     %   See matRad_DoseConstraint for interface description
0004     %
0005     % use log sum exp approximation, see appendix A in
0006     % http://scitation.aip.org/content/aapm/journal/medphys/41/8/10.1118/1.4883837
0007     %
0008     % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0009     %
0010     % Copyright 2020 the matRad development team.
0011     %
0012     % This file is part of the matRad project. It is subject to the license
0013     % terms in the LICENSE file found in the top-level directory of this
0014     % distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
0015     % of the matRad project, including this file, may be copied, modified,
0016     % propagated, or distributed except according to the terms contained in the
0017     % LICENSE file.
0018     %
0019     % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0020     
0021     properties (Constant)
0022         name = 'Min/Max dose constraint';
0023         parameterNames = {'d^{min}', 'd^{max}','method'};
0024         parameterTypes = {'dose','dose',{'approx','voxelwise'}};
0025     end
0026     
0027     properties
0028         parameters = {0,30,1};
0029         epsilon = 1e-3; %slack parameter for the logistic approximation
0030     end
0031     
0032     methods
0033         function obj = matRad_MinMaxDose(minDose,maxDose,method)
0034             
0035             %If we have a struct in first argument
0036             if nargin == 1 && isstruct(minDose)
0037                 inputStruct = minDose;
0038                 initFromStruct = true;
0039             else
0040                 initFromStruct = false;
0041                 inputStruct = [];
0042             end
0043             
0044             %Call Superclass Constructor (for struct initialization)
0045             obj@DoseConstraints.matRad_DoseConstraint(inputStruct);
0046             
0047             %now handle initialization from other parameters
0048             if ~initFromStruct
0049                 
0050                 if nargin < 3 || ~ischar(method)
0051                     method = 'approx';
0052                 end
0053                 
0054                 methodIx = find(strcmp(method,obj.parameterTypes{3}));
0055                 
0056                 if isempty(methodIx) || numel(methodIx) > 1
0057                     methodIx = 1;
0058                     msg = ['Dose Constraint method can only be ', strjoin(obj.parameterTypes{3},' or '), '! Using method ''', obj.parameterTypes{3}{methodIx}, '''.'];
0059                     warning(msg);
0060                 end
0061                 
0062                 obj.parameters{3} = methodIx;
0063                 
0064                 if nargin >= 2 && isscalar(maxDose)
0065                     obj.parameters{2} = maxDose;
0066                 end
0067                 
0068                 if nargin >= 1 && isscalar(minDose)
0069                     obj.parameters{1} = minDose;
0070                 end
0071             end
0072         end
0073         
0074         %Overloads the struct function to add constraint specific
0075         %parameters
0076         function s = struct(obj)
0077             s = struct@DoseConstraints.matRad_DoseConstraint(obj);
0078             s.epsilon = obj.epsilon;
0079         end
0080         
0081         function cu = upperBounds(obj,n)
0082             switch obj.parameters{3}
0083                 case 1 %logsumexp approx
0084                     %Validate parameters
0085                     if obj.parameters{1} <= 0 && isinf(obj.parameters{2}) %Constraint doesn't make sense (min = 0 & max = Inf)
0086                         cu = [];
0087                     elseif obj.parameters{2} == Inf %Only min dose
0088                         cu = Inf;
0089                     elseif obj.parameters{1} <= 0 %Only max dose
0090                         cu = obj.parameters{2};
0091                     else %both are set sensible
0092                         cu = [Inf; obj.parameters{2}];
0093                     end
0094                 case 2 %voxelwise
0095                     cu = obj.parameters{2}*ones(n,1);
0096                 otherwise
0097                     error(['Min/max dose constraint evaluation method not known!']);
0098             end
0099             %cu = [Inf; obj.parameters{2}];
0100         end
0101         function cl = lowerBounds(obj,n)
0102             switch obj.parameters{3}
0103                 case 1 %logsumexp approx
0104                     if obj.parameters{1} <= 0 && isinf(obj.parameters{2})
0105                         cl = [];
0106                     elseif obj.parameters{2} == Inf
0107                         cl = obj.parameters{1};
0108                     elseif obj.parameters{1} <= 0
0109                         cl = 0;
0110                     else
0111                         cl = [obj.parameters{1}; 0];
0112                     end
0113                 case 2
0114                     cl = obj.parameters{1}*ones(n,1);
0115                 otherwise
0116                     matRad_cfg = MatRad_Config.instance();
0117                     matRad_cfg.dispError('Min/max dose constraint evaluation method not known!');
0118             end
0119         end
0120         
0121         function jStruct = getDoseConstraintJacobianStructure(obj,n)
0122             switch obj.parameters{3}
0123                 case 1 %logsumexp approx
0124                     %Validate parameters
0125                     if obj.parameters{1} <= 0 && isinf(obj.parameters{2}) %Constraint doesn't make sense (min = 0 & max = Inf)
0126                         jStruct = ones(n,0);
0127                     elseif obj.parameters{1} > 0 && isfinite(obj.parameters{2}) %both are set sensible
0128                         jStruct = ones(n,2);
0129                     else %Only min or max dose
0130                         jStruct = ones(n,1);
0131                     end
0132                     %jStruct = ones(n,2);
0133                 case 2
0134                     jStruct = speye(n);
0135                 otherwise
0136                     matRad_cfg = MatRad_Config.instance();
0137                     matRad_cfg.dispError('Min/max dose constraint evaluation method not known!');
0138             end
0139             
0140         end
0141         
0142         %% Calculates the Constraint Function value
0143         function cDose = computeDoseConstraintFunction(obj,dose)
0144             %cDose(2) = dose_max + obj.epsilon * log( sum(exp((dose - dose_max)/obj.epsilon)));
0145             %cDose(1) = dose_min - obj.epsilon * log( sum(exp((dose_min - dose)/obj.epsilon)));
0146             switch obj.parameters{3}
0147                 case 1 %logsumexp approx
0148                     cDose = obj.computeDoseConstraintFunctionLogSumExp(dose);
0149                 case 2
0150                     cDose = obj.computeDoseConstraintFunctionVoxelwise(dose);
0151                 otherwise
0152                     matRad_cfg = MatRad_Config.instance();
0153                     matRad_cfg.dispError('Min/max dose constraint evaluation method not known!');
0154             end
0155         end
0156         
0157         %% Calculates the Constraint jacobian
0158         function cDoseJacob  = computeDoseConstraintJacobian(obj,dose)
0159             switch obj.parameters{3}
0160                 case 1 %logsumexp approx
0161                     cDoseJacob = obj.computeDoseConstraintJacobianLogSumExp(dose);
0162                 case 2
0163                     cDoseJacob = obj.computeDoseConstraintJacobianVoxelwise(dose);
0164                 otherwise
0165                     matRad_cfg = MatRad_Config.instance();
0166                     matRad_cfg.dispError('Min/max dose constraint evaluation method not known!');
0167             end
0168         end
0169     end
0170     
0171     methods (Access = private)
0172         % LogSumExp Approximation
0173         function cDose = computeDoseConstraintFunctionLogSumExp(obj,dose)
0174             dose_min = min(dose);
0175             dose_max = max(dose);
0176             
0177             %Validate parameters
0178             if obj.parameters{1} <= 0 && isinf(obj.parameters{2}) %Constraint doesn't make sense (min = 0 & max = Inf)
0179                 cDose = [];
0180             elseif obj.parameters{2} == Inf %Only min dose
0181                 cDose = dose_min - obj.epsilon * log( sum(exp((dose_min - dose)/obj.epsilon)));
0182             elseif obj.parameters{1} <= 0 %Only max dose
0183                 cDose = dose_max + obj.epsilon * log( sum(exp((dose - dose_max)/obj.epsilon)));
0184             else %both are set sensible
0185                 cDose(2,1) = dose_max + obj.epsilon * log( sum(exp((dose - dose_max)/obj.epsilon)));
0186                 cDose(1,1) = dose_min - obj.epsilon * log( sum(exp((dose_min - dose)/obj.epsilon)));
0187             end
0188             
0189         end
0190         function cDoseJacob  = computeDoseConstraintJacobianLogSumExp(obj,dose)
0191             %Validate parameters
0192             if obj.parameters{1} <= 0 && isinf(obj.parameters{2}) %Constraint doesn't make sense (min = 0 & max = Inf)
0193                 cDoseJacob = [];
0194             elseif obj.parameters{2} == Inf %Only min dose
0195                 cDoseJacob(:,1) = exp( (min(dose)-dose)/obj.epsilon );
0196                 cDoseJacob(:,1) = cDoseJacob(:,1)/sum(cDoseJacob(:,1));
0197             elseif obj.parameters{1} <= 0 %Only max dose
0198                 cDoseJacob(:,1) = exp( (dose-max(dose))/obj.epsilon );
0199                 cDoseJacob(:,1) = cDoseJacob(:,1)/sum(cDoseJacob(:,1));
0200             else %both are set sensible
0201                 cDoseJacob(:,1) = exp( (min(dose)-dose)/obj.epsilon );
0202                 cDoseJacob(:,1) = cDoseJacob(:,1)/sum(cDoseJacob(:,1));
0203                 
0204                 cDoseJacob(:,2) = exp( (dose-max(dose))/obj.epsilon );
0205                 cDoseJacob(:,2) = cDoseJacob(:,2)/sum(cDoseJacob(:,2));
0206             end
0207             
0208             
0209         end
0210         
0211         %Exact voxel-wise
0212         function cDose = computeDoseConstraintFunctionVoxelwise(obj,dose)
0213             cDose = dose;
0214         end
0215         function cDoseJacob  = computeDoseConstraintJacobianVoxelwise(obj,dose)
0216             cDoseJacob = speye(numel(dose),numel(dose));
0217         end
0218     end
0219     
0220 end
0221 
0222