From 8e78438f06572543e09968d293288f5e14feec58 Mon Sep 17 00:00:00 2001 From: Niklas Wahl Date: Thu, 26 Mar 2026 18:22:45 +0100 Subject: [PATCH 1/2] fixes the GriddedScenarios implementation when the number of grid points is 1, which will reduce to the nominal scenario in this dimension. Also provides corresponding tests for this case --- .../matRad_GriddedScenariosAbstract.m | 305 +++++++-------- matRad/scenarios/matRad_WorstCaseScenarios.m | 103 ++++-- test/scenarios/test_wcScenarios.m | 350 ++++++++++-------- 3 files changed, 397 insertions(+), 361 deletions(-) diff --git a/matRad/scenarios/matRad_GriddedScenariosAbstract.m b/matRad/scenarios/matRad_GriddedScenariosAbstract.m index 8f34f3f9d..339041b68 100644 --- a/matRad/scenarios/matRad_GriddedScenariosAbstract.m +++ b/matRad/scenarios/matRad_GriddedScenariosAbstract.m @@ -1,58 +1,61 @@ classdef (Abstract) matRad_GriddedScenariosAbstract < matRad_ScenarioModel - %UNTITLED Summary of this class goes here + % UNTITLED Summary of this class goes here % Detailed explanation goes here properties (AbortSet = true) - %includeNominalScenario = true; - combinations = 'none'; %Can be 'none', 'shift', 'all' to ontrol creation of worst case combinations - combineRange = true; %Wether to treat absolute & relative range as one shift or as separate scenarios + % includeNominalScenario = true; + combinations = 'none' % Can be 'none', 'shift', 'all' to control creation of worst case combinations + combineRange = true % Whether to treat absolute & relative range as one shift or as separate scenarios end - - %Each subclass needs to define how many gridpoints it uses and if this - %can be set or not + + % Each subclass needs to define how many gridpoints it uses and if this + % can be set or not properties (Abstract) - numOfSetupGridPoints; - numOfRangeGridPoints; + numOfSetupGridPoints + numOfRangeGridPoints end properties (Constant) - validCombinationTypes = {'all','none','shift'}; + validCombinationTypes = {'all', 'none', 'shift'} end methods - function this = matRad_GriddedScenariosAbstract(ct) - - if nargin == 0 + + function this = matRad_GriddedScenariosAbstract(ct) + + if nargin == 0 superclassArgs = {}; else superclassArgs = {ct}; end - + this@matRad_ScenarioModel(superclassArgs{:}); - %TODO: We could do this automatically in the superclass - %Octave 5 has a bug there and throws an error - %this.updateScenarios(); + % TODO: We could do this automatically in the superclass + % Octave 5 has a bug there and throws an error + % this.updateScenarios(); end - + %% set methods - function set.combineRange(this,combineRange_) - valid = isscalar(combineRange_) && (isnumeric(combineRange_) || islogical(combineRange_)); - if ~valid + function set.combineRange(this, combineRange) + valid = isscalar(combineRange) && (isnumeric(combineRange) || islogical(combineRange)); + if ~valid matRad_cfg = MatRad_Config.instance(); matRad_cfg.dispError('Invalid value for combineRange! Needs to be a boolean / logical value!'); end - this.combineRange = combineRange_; + this.combineRange = combineRange; this.updateScenarios(); end - function set.combinations(this,combinations_) - valid = any(strcmp(combinations_,this.validCombinationTypes)); - if ~valid + function set.combinations(this, combinations) + valid = any(strcmp(combinations, this.validCombinationTypes)); + if ~valid matRad_cfg = MatRad_Config.instance(); - matRad_cfg.dispError('Invalid value for combinations! Needs to be one of the strings %s!',strjoin(this.validCombinationTypes,' / ')); + matRad_cfg.dispError( ... + 'Invalid value for combinations! Needs to be one of the strings %s!', ... + strjoin(this.validCombinationTypes, ' / ')); end - this.combinations = combinations_; + this.combinations = combinations; this.updateScenarios(); end @@ -60,211 +63,191 @@ matRad_cfg = MatRad_Config.instance(); % - this.numOfCtScen = size(this.ctScenProb,1); + this.numOfCtScen = size(this.ctScenProb, 1); - %Get the maximum, i.e., worst case shifts + % Get the maximum, i.e., worst case shifts wcSetupShifts = this.wcSigma * this.shiftSD; - + %% Create gridded setup shifts - %Create grid vectors for setup shifts - setupShiftGrid = zeros(this.numOfSetupGridPoints,numel(wcSetupShifts)); - %{ - if mod(this.numOfSetupGridPoints,2) == 0 && this.includeNominalScenario - matRad_cfg.dispWarning('Obtaining Setup Shifts: Including the nominal scenario with even number of grid points creates asymmetrical shifts!'); - end - %} + % Create grid vectors for setup shifts + setupShiftGrid = zeros(this.numOfSetupGridPoints, numel(wcSetupShifts)); for i = 1:numel(wcSetupShifts) - setupShiftGrid(:,i) = linspace(-wcSetupShifts(i),wcSetupShifts(i),this.numOfSetupGridPoints); - %{ - if this.includeNominalScenario - - [~,ix] = min(abs(setupShiftGrid(:,i))); - setupShiftGrid(ix,i) = 0; - end - %} + if size(setupShiftGrid, 1) > 1 + setupShiftGrid(:, i) = linspace(-wcSetupShifts(i), wcSetupShifts(i), this.numOfSetupGridPoints); + else + setupShiftGrid(1, i) = 0; + end end - - %Now create vector of all shifts for different combinatorial - %settings + + % Now create vector of all shifts for different combinatorial + % settings switch this.combinations case 'none' - %Create independent shifts + % Create independent shifts griddedSetupShifts = []; - for i=1:size(setupShiftGrid,2) - tmpGrid = zeros(size(setupShiftGrid,1),3); - tmpGrid(:,i) = setupShiftGrid(:,i); + for i = 1:size(setupShiftGrid, 2) + tmpGrid = zeros(size(setupShiftGrid, 1), 3); + tmpGrid(:, i) = setupShiftGrid(:, i); griddedSetupShifts = [griddedSetupShifts; tmpGrid]; - end - case {'shift','all'} - [X,Y,Z] = meshgrid(setupShiftGrid(:,1),setupShiftGrid(:,2),setupShiftGrid(:,3)); - griddedSetupShifts = [X(:), Y(:), Z(:)]; + end + case {'shift', 'all'} + [X, Y, Z] = meshgrid(setupShiftGrid(:, 1), setupShiftGrid(:, 2), setupShiftGrid(:, 3)); + griddedSetupShifts = [X(:), Y(:), Z(:)]; otherwise matRad_cfg.dispError('Invalid value for combinations! This sanity check should never be reached!'); end griddedSetupShifts = matRad_ImportanceScenarios.uniqueStableRowsCompat(griddedSetupShifts); - shiftNomScenIx = find(all(griddedSetupShifts == zeros(1,3),2)); - + shiftNomScenIx = find(all(griddedSetupShifts == zeros(1, 3), 2)); + if ~isempty(shiftNomScenIx) %|| this.includeNominalScenario if ~isempty(shiftNomScenIx) - griddedSetupShifts(shiftNomScenIx,:) = []; + griddedSetupShifts(shiftNomScenIx, :) = []; end griddedSetupShifts = [0 0 0; griddedSetupShifts]; end - - this.totNumShiftScen = size(griddedSetupShifts,1); - + + this.totNumShiftScen = size(griddedSetupShifts, 1); + %% Create gridded range shifts - %Obtain worst case range shifts - wcRangeShifts = this.wcSigma * [this.rangeAbsSD this.rangeRelSD./100]; - - rangeShiftGrid = zeros(this.numOfRangeGridPoints,numel(wcRangeShifts)); - %{ - if mod(this.numOfRangeGridPoints,2) == 0 && this.includeNominalScenario - matRad_cfg.dispWarning('Obtaining Range Shifts: Including the nominal scenario with even number of grid points creates asymmetrical shifts!'); - end - %} + % Obtain worst case range shifts + wcRangeShifts = this.wcSigma * [this.rangeAbsSD this.rangeRelSD ./ 100]; + + rangeShiftGrid = zeros(this.numOfRangeGridPoints, numel(wcRangeShifts)); for i = 1:numel(wcRangeShifts) - rangeShiftGrid(:,i) = linspace(-wcRangeShifts(i),wcRangeShifts(i),this.numOfRangeGridPoints); - - %{ - if this.includeNominalScenario - [~,ix] = min(abs(rangeShiftGrid(:,i))); - rangeShiftGrid(ix,i) = 0; + if size(rangeShiftGrid, 1) > 1 + rangeShiftGrid(:, i) = linspace(-wcRangeShifts(i), wcRangeShifts(i), this.numOfRangeGridPoints); + else + rangeShiftGrid(1, i) = 0; end - %} end if this.combineRange griddedRangeShifts = rangeShiftGrid; - else - [rngAbs,rngRel] = meshgrid(rangeShiftGrid(:,1),rangeShiftGrid(:,2)); - griddedRangeShifts = [rngAbs(:),rngRel(:)]; + else + [rngAbs, rngRel] = meshgrid(rangeShiftGrid(:, 1), rangeShiftGrid(:, 2)); + griddedRangeShifts = [rngAbs(:), rngRel(:)]; end - %Remove duplicate scenarios and update number of shifts - griddedRangeShifts = this.uniqueStableRowsCompat(griddedRangeShifts); + % Remove duplicate scenarios and update number of shifts + griddedRangeShifts = this.uniqueStableRowsCompat(griddedRangeShifts); + + rangeNomScenIx = find(all(griddedRangeShifts == zeros(1, 2), 2)); - rangeNomScenIx = find(all(griddedRangeShifts == zeros(1,2),2)); - if ~isempty(rangeNomScenIx) %|| this.includeNominalScenario if ~isempty(rangeNomScenIx) - griddedRangeShifts(rangeNomScenIx,:) = []; + griddedRangeShifts(rangeNomScenIx, :) = []; end griddedRangeShifts = [0 0; griddedRangeShifts]; end - this.totNumRangeScen = size(griddedRangeShifts,1); - - %Aggregate scenarios + this.totNumRangeScen = size(griddedRangeShifts, 1); + + % Aggregate scenarios switch this.combinations - case {'none','shift'} - scenarios = zeros(this.totNumShiftScen + this.totNumRangeScen,5); - scenarios(1:this.totNumShiftScen,1:3) = griddedSetupShifts; - scenarios(this.totNumShiftScen+1:this.totNumShiftScen+this.totNumRangeScen,4:5) = griddedRangeShifts; + case {'none', 'shift'} + scenarios = zeros(this.totNumShiftScen + this.totNumRangeScen, 5); + scenarios(1:this.totNumShiftScen, 1:3) = griddedSetupShifts; + scenarios(this.totNumShiftScen + 1:this.totNumShiftScen + this.totNumRangeScen, 4:5) = griddedRangeShifts; - %create the linear mask of scenarios - linearMaskTmp = ones(size(scenarios,1),3); - linearMaskTmp(1:this.totNumShiftScen,2) = (1:this.totNumShiftScen)'; - linearMaskTmp(this.totNumShiftScen+1:this.totNumShiftScen+this.totNumRangeScen,3) = (1:this.totNumRangeScen)'; + % create the linear mask of scenarios + linearMaskTmp = ones(size(scenarios, 1), 3); + linearMaskTmp(1:this.totNumShiftScen, 2) = (1:this.totNumShiftScen)'; + linearMaskTmp(this.totNumShiftScen + 1:this.totNumShiftScen + this.totNumRangeScen, 3) = (1:this.totNumRangeScen)'; - [scenarios,ia] = matRad_ImportanceScenarios.uniqueStableRowsCompat(scenarios); - linearMaskTmp = linearMaskTmp(ia,:); + [scenarios, ia] = matRad_ImportanceScenarios.uniqueStableRowsCompat(scenarios); + linearMaskTmp = linearMaskTmp(ia, :); case 'all' - %Prepare scenario matrix by replicating shifts - %with the number of range scenarios - scenarios = repmat(griddedSetupShifts,this.totNumRangeScen,1); - scenarios = [scenarios zeros(size(scenarios,1),2)]; - - %create the linear mask of scenarios - linearMaskTmp = ones(size(scenarios,1),3); - for r = 1:this.totNumRangeScen - offset = (r-1)*this.totNumShiftScen; + % Prepare scenario matrix by replicating shifts + % with the number of range scenarios + scenarios = repmat(griddedSetupShifts, this.totNumRangeScen, 1); + scenarios = [scenarios zeros(size(scenarios, 1), 2)]; + + % create the linear mask of scenarios + linearMaskTmp = ones(size(scenarios, 1), 3); + for r = 1:this.totNumRangeScen + offset = (r - 1) * this.totNumShiftScen; ixR = (offset + 1):(offset + this.totNumShiftScen); - scenarios(ixR,4:5) = repmat(griddedRangeShifts(r,:),this.totNumShiftScen,1); + scenarios(ixR, 4:5) = repmat(griddedRangeShifts(r, :), this.totNumShiftScen, 1); - %Set linear mask - linearMaskTmp(ixR,2) = (1:this.totNumShiftScen)'; - linearMaskTmp(ixR,3) = r; + % Set linear mask + linearMaskTmp(ixR, 2) = (1:this.totNumShiftScen)'; + linearMaskTmp(ixR, 3) = r; end - - %create the linear mask of scenarios - [scenarios,ia] = matRad_ImportanceScenarios.uniqueStableRowsCompat(scenarios); - linearMaskTmp = linearMaskTmp(ia,:); + % create the linear mask of scenarios + [scenarios, ia] = matRad_ImportanceScenarios.uniqueStableRowsCompat(scenarios); + linearMaskTmp = linearMaskTmp(ia, :); otherwise matRad_cfg.dispError('Invalid value for combinations! This sanity check should never be reached!'); end - %if ~this.includeNominalScenario - % nomScen = all(scenarios == zeros(1,5),2); - % scenarios(nomScen,:) = []; - % linearMaskTmp(nomScen,:) = []; - %end - - %Handle 4D phases - phases = repmat(this.ctScenProb(:,1)',size(scenarios,1),1); + % Handle 4D phases + phases = repmat(this.ctScenProb(:, 1)', size(scenarios, 1), 1); phases = phases(:); - scenarios = horzcat(phases, repmat(scenarios,[this.numOfCtScen 1])); - linearMaskTmp = repmat(linearMaskTmp,this.numOfCtScen,1); - linearMaskTmp(:,1) = phases; + scenarios = horzcat(phases, repmat(scenarios, [this.numOfCtScen 1])); + linearMaskTmp = repmat(linearMaskTmp, this.numOfCtScen, 1); + linearMaskTmp(:, 1) = phases; this.ctScenIx = phases; - %Finalize meta information - this.totNumScen = size(scenarios,1); + % Finalize meta information + this.totNumScen = size(scenarios, 1); + + this.relRangeShift = scenarios(:, 6); + this.absRangeShift = scenarios(:, 5); + this.isoShift = scenarios(:, 2:4); - this.relRangeShift = scenarios(:,6); - this.absRangeShift = scenarios(:,5); - this.isoShift = scenarios(:,2:4); - this.maxAbsRangeShift = max(this.absRangeShift); this.maxRelRangeShift = max(this.relRangeShift); - this.scenMask = false(this.numOfAvailableCtScen,this.totNumShiftScen,this.totNumRangeScen); - + this.scenMask = false(this.numOfAvailableCtScen, this.totNumShiftScen, this.totNumRangeScen); + this.scenForProb = scenarios; this.linearMask = linearMaskTmp; - - maskIx = sub2ind(size(this.scenMask),linearMaskTmp(:,1),linearMaskTmp(:,2),linearMaskTmp(:,3)); + + maskIx = sub2ind(size(this.scenMask), linearMaskTmp(:, 1), linearMaskTmp(:, 2), linearMaskTmp(:, 3)); this.scenMask(maskIx) = true; - %Get Scenario probability - %First, we use the Gaussian Uncertainty model for range and - %setup - Sigma = diag([this.shiftSD,this.rangeAbsSD,this.rangeRelSD./100].^2); - d = size(Sigma,1); - [cs,p] = chol(Sigma); - - tmpScenProb = (2*pi)^(-d/2) * exp(-0.5*sum((scenarios(:,2:end)/cs).^2, 2)) / prod(diag(cs)); - - %Now we combine with the 4D ct phase probabilities (multiply) - tmpPhaseProb = arrayfun(@(phase) this.ctScenProb(find(this.ctScenProb(:,1) == phase),2),phases); - - %Finalize probabilities + % Get Scenario probability + % First, we use the Gaussian Uncertainty model for range and + % setup + Sigma = diag([this.shiftSD, this.rangeAbsSD, this.rangeRelSD ./ 100].^2); + d = size(Sigma, 1); + [cs, p] = chol(Sigma); + + tmpScenProb = (2 * pi)^(-d / 2) * exp(-0.5 * sum((scenarios(:, 2:end) / cs).^2, 2)) / prod(diag(cs)); + + % Now we combine with the 4D ct phase probabilities (multiply) + tmpPhaseProb = arrayfun(@(phase) this.ctScenProb(find(this.ctScenProb(:, 1) == phase), 2), phases); + + % Finalize probabilities this.scenProb = tmpPhaseProb .* tmpScenProb; - this.scenWeight = this.scenProb./sum(this.scenProb); + this.scenWeight = this.scenProb ./ sum(this.scenProb); - %TODO: Discard scenarios with probability 0? + % TODO: Discard scenarios with probability 0? end + end methods (Static) - function [uniqueStableRows,ia] = uniqueStableRowsCompat(values) - %This is a compatability wrapper to call unique without sorting - + + function [uniqueStableRows, ia] = uniqueStableRowsCompat(values) + % This is a compatibility wrapper to call unique without sorting + matRad_cfg = MatRad_Config.instance(); - + if matRad_cfg.isOctave - %https://stackoverflow.com/questions/37828894/ - [~,ia,~] = unique(values,'rows','first'); + % https://stackoverflow.com/questions/37828894/ + [~, ia, ~] = unique(values, 'rows', 'first'); ia = sort(ia); - uniqueStableRows = values(ia,:); + uniqueStableRows = values(ia, :); else - [uniqueStableRows,ia] = unique(values,'rows','stable'); + [uniqueStableRows, ia] = unique(values, 'rows', 'stable'); end end + end -end \ No newline at end of file +end diff --git a/matRad/scenarios/matRad_WorstCaseScenarios.m b/matRad/scenarios/matRad_WorstCaseScenarios.m index 434cdf6c4..084d9288d 100644 --- a/matRad/scenarios/matRad_WorstCaseScenarios.m +++ b/matRad/scenarios/matRad_WorstCaseScenarios.m @@ -1,60 +1,83 @@ classdef matRad_WorstCaseScenarios < matRad_GriddedScenariosAbstract -% matRad_WorstCaseScenarios -% Implements single worst-case shifts per dimension.% -% -% constructor -% matRad_WorstCaseScenarios() -% matRad_WorstCaseScenarios(ct) -% -% input: -% ct: ct cube -% -% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% -% Copyright 2022-2026 the matRad development team. -% -% This file is part of the matRad project. It is subject to the license -% terms in the LICENSE file found in the top-level directory of this -% distribution and at https://github.com/e0404/matRad/LICENSE.md. No part -% of the matRad project, including this file, may be copied, modified, -% propagated, or distributed except according to the terms contained in the -% LICENSE file. -% -% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - - properties (SetAccess=protected) - shortName = 'wcScen'; - name = 'Worst Case Scenarios'; + % matRad_WorstCaseScenarios + % Implements single worst-case shifts per dimension.% + % + % constructor + % matRad_WorstCaseScenarios() + % matRad_WorstCaseScenarios(ct) + % + % input: + % ct: ct cube + % + % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % + % Copyright 2022-2026 the matRad development team. + % + % This file is part of the matRad project. It is subject to the license + % terms in the LICENSE file found in the top-level directory of this + % distribution and at https://github.com/e0404/matRad/LICENSE.md. No part + % of the matRad project, including this file, may be copied, modified, + % propagated, or distributed except according to the terms contained in the + % LICENSE file. + % + % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + properties (SetAccess = protected) + shortName = 'wcScen' + name = 'Worst Case Scenarios' end properties (Hidden) - numOfSetupGridPoints = 3; - numOfRangeGridPoints = 3; + numOfSetupGridPoints = 3 + numOfRangeGridPoints = 3 end - + methods - function this = matRad_WorstCaseScenarios(ct) - if nargin == 0 + + function this = matRad_WorstCaseScenarios(ct) + if nargin == 0 superclassArgs = {}; else superclassArgs = {ct}; end - + this@matRad_GriddedScenariosAbstract(superclassArgs{:}); - %TODO: We could do this automatically in the superclass - %Octave 5 has a bug there and throws an error + % TODO: We could do this automatically in the superclass + % Octave 5 has a bug there and throws an error this.updateScenarios(); end - + function scenarios = updateScenarios(this) - %Use the static gridded shift function from - %ImportanceScenarios. We set inclusion of nominal scenarios to - %false and handle it automatically via the grid point number - scenarios = this.updateScenarios@matRad_GriddedScenariosAbstract(); + % Use the static gridded shift function from + % ImportanceScenarios. We set inclusion of nominal scenarios to + % false and handle it automatically via the grid point number + scenarios = this.updateScenarios@matRad_GriddedScenariosAbstract(); end + end + methods + + function set.numOfSetupGridPoints(this, numOfSetupGridPoints) + if numOfSetupGridPoints ~= 3 && numOfSetupGridPoints ~= 1 + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispError([mfilename('class') ... + ' only supports 1 (only nominal grid point) or 3 grid points (nominal + worst case) for setup uncertainties!']); + end + this.numOfSetupGridPoints = numOfSetupGridPoints; + this.updateScenarios(); + end + + function set.numOfRangeGridPoints(this, numOfRangeGridPoints) + if numOfRangeGridPoints ~= 3 && numOfRangeGridPoints ~= 1 + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispError([mfilename('class') ... + ' only supports 1 (only nominal grid point) or 3 grid points (nominal + worst case) for range uncertainties!']); + end + this.numOfRangeGridPoints = numOfRangeGridPoints; + this.updateScenarios(); + end + end end diff --git a/test/scenarios/test_wcScenarios.m b/test/scenarios/test_wcScenarios.m index def02eb5a..3d646db99 100644 --- a/test/scenarios/test_wcScenarios.m +++ b/test/scenarios/test_wcScenarios.m @@ -1,178 +1,208 @@ function test_suite = test_wcScenarios -test_functions=localfunctions(); +test_functions = localfunctions(); initTestSuite; function test_worstCaseScenarioConstructor - scenario = matRad_WorstCaseScenarios(); - assertTrue(isa(scenario, 'matRad_WorstCaseScenarios')); - assertTrue(isa(scenario, 'matRad_GriddedScenariosAbstract')); - assertTrue(isa(scenario, 'matRad_ScenarioModel')); - assertEqual(scenario.shortName, 'wcScen'); - %Test correct standard values & sizes - assertEqual(scenario.ctScenProb, [1 1]); - assertEqual(scenario.numOfCtScen, 1); - assertEqual(scenario.totNumScen, 9); - assertEqual(scenario.totNumShiftScen, 7); - assertEqual(scenario.totNumRangeScen, 3); - assertEqual(size(scenario.relRangeShift),[scenario.totNumScen,1]); - assertEqual(size(scenario.absRangeShift),[scenario.totNumScen,1]); - assertEqual(size(scenario.isoShift),[scenario.totNumScen,3]); - assertEqual(scenario.relRangeShift, [zeros(7,1); -0.035; 0.035]); - assertEqual(scenario.absRangeShift, [zeros(7,1); -1; 1]); - %assertEqual(scenario.isoShift, 2.25 * ones(1,3)); - assertEqual(scenario.maxAbsRangeShift, max(scenario.absRangeShift)); - assertEqual(scenario.maxRelRangeShift, max(scenario.relRangeShift)); - assertEqual(size(scenario.scenMask), [scenario.numOfCtScen,scenario.totNumShiftScen,scenario.totNumRangeScen]); - %assertEqual(scenario.scenMask, true(1,1,1)); - assertEqual(size(scenario.linearMask), [scenario.totNumScen,3]); - - [tmp(:,1),tmp(:,2),tmp(:,3)] = ind2sub(size(scenario.scenMask),find(scenario.scenMask)); - assertEqual(tmp,scenario.linearMask); - - %assertEqual(ind2sub(find())) - %assertEqual(scenario.linearMask, [1 1 1]); - assertElementsAlmostEqual(scenario.scenProb,helper_mvarGauss(scenario)); - - tmp = [scenario.ctScenIx scenario.isoShift scenario.absRangeShift scenario.relRangeShift]; - assertEqual(scenario.scenForProb,tmp); - assertEqual(scenario.scenWeight,scenario.scenProb./sum(scenario.scenProb)); +scenario = matRad_WorstCaseScenarios(); +assertTrue(isa(scenario, 'matRad_WorstCaseScenarios')); +assertTrue(isa(scenario, 'matRad_GriddedScenariosAbstract')); +assertTrue(isa(scenario, 'matRad_ScenarioModel')); +assertEqual(scenario.shortName, 'wcScen'); +% Test correct standard values & sizes +assertEqual(scenario.ctScenProb, [1 1]); +assertEqual(scenario.numOfCtScen, 1); +assertEqual(scenario.totNumScen, 9); +assertEqual(scenario.totNumShiftScen, 7); +assertEqual(scenario.totNumRangeScen, 3); +assertEqual(size(scenario.relRangeShift), [scenario.totNumScen, 1]); +assertEqual(size(scenario.absRangeShift), [scenario.totNumScen, 1]); +assertEqual(size(scenario.isoShift), [scenario.totNumScen, 3]); +assertEqual(scenario.relRangeShift, [zeros(7, 1); -0.035; 0.035]); +assertEqual(scenario.absRangeShift, [zeros(7, 1); -1; 1]); +% assertEqual(scenario.isoShift, 2.25 * ones(1,3)); +assertEqual(scenario.maxAbsRangeShift, max(scenario.absRangeShift)); +assertEqual(scenario.maxRelRangeShift, max(scenario.relRangeShift)); +assertEqual(size(scenario.scenMask), [scenario.numOfCtScen, scenario.totNumShiftScen, scenario.totNumRangeScen]); +% assertEqual(scenario.scenMask, true(1,1,1)); +assertEqual(size(scenario.linearMask), [scenario.totNumScen, 3]); + +[tmp(:, 1), tmp(:, 2), tmp(:, 3)] = ind2sub(size(scenario.scenMask), find(scenario.scenMask)); +assertEqual(tmp, scenario.linearMask); + +% assertEqual(ind2sub(find())) +% assertEqual(scenario.linearMask, [1 1 1]); +assertElementsAlmostEqual(scenario.scenProb, helper_mvarGauss(scenario)); + +tmp = [scenario.ctScenIx scenario.isoShift scenario.absRangeShift scenario.relRangeShift]; +assertEqual(scenario.scenForProb, tmp); +assertEqual(scenario.scenWeight, scenario.scenProb ./ sum(scenario.scenProb)); function test_worstCaseScenarioConstructorWithCt - n = 5; - ct = struct('numOfCtScen',n); - - scenario = matRad_WorstCaseScenarios(ct); - assertTrue(isa(scenario, 'matRad_WorstCaseScenarios')); - assertTrue(isa(scenario, 'matRad_GriddedScenariosAbstract')); - assertTrue(isa(scenario, 'matRad_ScenarioModel')); - assertEqual(scenario.shortName, 'wcScen'); - %Test correct standard values & sizes - assertEqual(scenario.ctScenProb, [(1:n)' ones(n,1)./n]); - assertEqual(scenario.numOfCtScen, n); - assertEqual(scenario.totNumScen, 9*n); - assertEqual(scenario.totNumShiftScen, 7); - assertEqual(scenario.totNumRangeScen, 3); - assertEqual(size(scenario.relRangeShift),[scenario.totNumScen,1]); - assertEqual(size(scenario.absRangeShift),[scenario.totNumScen,1]); - assertEqual(size(scenario.isoShift),[scenario.totNumScen,3]); - assertEqual(scenario.relRangeShift, repmat([zeros(7,1); -0.035; 0.035],n,1)); - assertEqual(scenario.absRangeShift, repmat([zeros(7,1); -1; 1],n,1)); - %assertEqual(scenario.isoShift, 2.25 * ones(1,3)); - assertEqual(scenario.maxAbsRangeShift, max(scenario.absRangeShift)); - assertEqual(scenario.maxRelRangeShift, max(scenario.relRangeShift)); - assertEqual(size(scenario.scenMask), [scenario.numOfCtScen,scenario.totNumShiftScen,scenario.totNumRangeScen]); - %assertEqual(scenario.scenMask, true(1,1,1)); - assertEqual(size(scenario.linearMask), [scenario.totNumScen,3]); - - tmpScenMask = permute(scenario.scenMask,[2 3 1]); - - [tmp(:,2),tmp(:,3),tmp(:,1)] = ind2sub(size(tmpScenMask),find(tmpScenMask)); - assertEqual(tmp,scenario.linearMask); - - %assertEqual(ind2sub(find())) - %assertEqual(scenario.linearMask, [1 1 1]); - assertElementsAlmostEqual(scenario.scenProb,helper_mvarGauss(scenario)); - - tmp = [scenario.ctScenIx scenario.isoShift scenario.absRangeShift scenario.relRangeShift]; - assertEqual(scenario.scenForProb,tmp); - assertEqual(scenario.scenWeight,scenario.scenProb./sum(scenario.scenProb)); - - -function test_worstCaseScenarioExtractSingleScenario - refScen = matRad_WorstCaseScenarios(); - - for scenNum = 1:refScen.totNumScen - scenario = refScen.extractSingleScenario(scenNum); - assertTrue(isa(scenario, 'matRad_NominalScenario')); - ctScenIx = refScen.ctScenIx(scenNum); - ctScenNum = find(ctScenIx == refScen.ctScenProb(:,1)); - assertEqual(scenario.ctScenProb, refScen.ctScenProb(ctScenNum,:)); - assertEqual(scenario.numOfCtScen, 1); - assertEqual(scenario.totNumScen, 1); - assertEqual(scenario.totNumShiftScen, 1); - assertEqual(scenario.totNumRangeScen, 1); - assertEqual(scenario.relRangeShift, refScen.relRangeShift(scenNum)); - assertEqual(scenario.absRangeShift, refScen.absRangeShift(scenNum)); - assertEqual(scenario.isoShift, refScen.isoShift(scenNum,:)); - assertEqual(scenario.maxAbsRangeShift, max(abs(refScen.absRangeShift(scenNum)))); - assertEqual(scenario.maxRelRangeShift, max(abs(refScen.relRangeShift(scenNum)))); - assertTrue(scenario.scenMask(ctScenIx,1,1)); - assertTrue(numel(find(scenario.scenMask)) == 1); - assertEqual(scenario.linearMask, [ctScenIx 1 1]); - assertElementsAlmostEqual(scenario.scenProb,helper_mvarGauss(scenario)); - assertEqual(scenario.scenForProb,refScen.scenForProb(scenNum,:)); - assertEqual(scenario.scenWeight, refScen.scenWeight(scenNum)); - end +n = 5; +ct = struct('numOfCtScen', n); + +scenario = matRad_WorstCaseScenarios(ct); +assertTrue(isa(scenario, 'matRad_WorstCaseScenarios')); +assertTrue(isa(scenario, 'matRad_GriddedScenariosAbstract')); +assertTrue(isa(scenario, 'matRad_ScenarioModel')); +assertEqual(scenario.shortName, 'wcScen'); +% Test correct standard values & sizes +assertEqual(scenario.ctScenProb, [(1:n)' ones(n, 1) ./ n]); +assertEqual(scenario.numOfCtScen, n); +assertEqual(scenario.totNumScen, 9 * n); +assertEqual(scenario.totNumShiftScen, 7); +assertEqual(scenario.totNumRangeScen, 3); +assertEqual(size(scenario.relRangeShift), [scenario.totNumScen, 1]); +assertEqual(size(scenario.absRangeShift), [scenario.totNumScen, 1]); +assertEqual(size(scenario.isoShift), [scenario.totNumScen, 3]); +assertEqual(scenario.relRangeShift, repmat([zeros(7, 1); -0.035; 0.035], n, 1)); +assertEqual(scenario.absRangeShift, repmat([zeros(7, 1); -1; 1], n, 1)); +% assertEqual(scenario.isoShift, 2.25 * ones(1,3)); +assertEqual(scenario.maxAbsRangeShift, max(scenario.absRangeShift)); +assertEqual(scenario.maxRelRangeShift, max(scenario.relRangeShift)); +assertEqual(size(scenario.scenMask), [scenario.numOfCtScen, scenario.totNumShiftScen, scenario.totNumRangeScen]); +% assertEqual(scenario.scenMask, true(1,1,1)); +assertEqual(size(scenario.linearMask), [scenario.totNumScen, 3]); + +tmpScenMask = permute(scenario.scenMask, [2 3 1]); + +[tmp(:, 2), tmp(:, 3), tmp(:, 1)] = ind2sub(size(tmpScenMask), find(tmpScenMask)); +assertEqual(tmp, scenario.linearMask); + +% assertEqual(ind2sub(find())) +% assertEqual(scenario.linearMask, [1 1 1]); +assertElementsAlmostEqual(scenario.scenProb, helper_mvarGauss(scenario)); + +tmp = [scenario.ctScenIx scenario.isoShift scenario.absRangeShift scenario.relRangeShift]; +assertEqual(scenario.scenForProb, tmp); +assertEqual(scenario.scenWeight, scenario.scenProb ./ sum(scenario.scenProb)); +function test_worstCaseScenarioExtractSingleScenario +refScen = matRad_WorstCaseScenarios(); + +for scenNum = 1:refScen.totNumScen + scenario = refScen.extractSingleScenario(scenNum); + assertTrue(isa(scenario, 'matRad_NominalScenario')); + ctScenIx = refScen.ctScenIx(scenNum); + ctScenNum = find(ctScenIx == refScen.ctScenProb(:, 1)); + assertEqual(scenario.ctScenProb, refScen.ctScenProb(ctScenNum, :)); + assertEqual(scenario.numOfCtScen, 1); + assertEqual(scenario.totNumScen, 1); + assertEqual(scenario.totNumShiftScen, 1); + assertEqual(scenario.totNumRangeScen, 1); + assertEqual(scenario.relRangeShift, refScen.relRangeShift(scenNum)); + assertEqual(scenario.absRangeShift, refScen.absRangeShift(scenNum)); + assertEqual(scenario.isoShift, refScen.isoShift(scenNum, :)); + assertEqual(scenario.maxAbsRangeShift, max(abs(refScen.absRangeShift(scenNum)))); + assertEqual(scenario.maxRelRangeShift, max(abs(refScen.relRangeShift(scenNum)))); + assertTrue(scenario.scenMask(ctScenIx, 1, 1)); + assertTrue(numel(find(scenario.scenMask)) == 1); + assertEqual(scenario.linearMask, [ctScenIx 1 1]); + assertElementsAlmostEqual(scenario.scenProb, helper_mvarGauss(scenario)); + assertEqual(scenario.scenForProb, refScen.scenForProb(scenNum, :)); + assertEqual(scenario.scenWeight, refScen.scenWeight(scenNum)); +end function test_worstCaseScenarioExtractSingleScenarioWithCtScen - n = 5; - scenNum = 1; - ct = struct('numOfCtScen',n); - refScen = matRad_WorstCaseScenarios(ct); - for scenNum = 1:refScen.totNumScen - scenario = refScen.extractSingleScenario(scenNum); - assertTrue(isa(scenario, 'matRad_NominalScenario')); - ctScenIx = refScen.ctScenIx(scenNum); - ctScenNum = find(ctScenIx == refScen.ctScenProb(:,1)); - assertEqual(scenario.ctScenProb, refScen.ctScenProb(ctScenNum,:)); - assertEqual(scenario.numOfCtScen, 1); - assertEqual(scenario.totNumScen, 1); - assertEqual(scenario.totNumShiftScen, 1); - assertEqual(scenario.totNumRangeScen, 1); - assertEqual(scenario.relRangeShift, refScen.relRangeShift(scenNum)); - assertEqual(scenario.absRangeShift, refScen.absRangeShift(scenNum)); - assertEqual(scenario.isoShift, refScen.isoShift(scenNum,:)); - assertEqual(scenario.maxAbsRangeShift, max(abs(refScen.absRangeShift(scenNum)))); - assertEqual(scenario.maxRelRangeShift, max(abs(refScen.relRangeShift(scenNum)))); - assertTrue(scenario.scenMask(ctScenIx,1,1)); - assertTrue(numel(find(scenario.scenMask)) == 1); - assertEqual(scenario.linearMask, [ctScenIx 1 1]); - assertElementsAlmostEqual(scenario.scenProb,helper_mvarGauss(scenario)); - assertEqual(scenario.scenForProb,refScen.scenForProb(scenNum,:)); - assertEqual(scenario.scenWeight, refScen.scenWeight(scenNum)); - end - +n = 5; +scenNum = 1; +ct = struct('numOfCtScen', n); +refScen = matRad_WorstCaseScenarios(ct); +for scenNum = 1:refScen.totNumScen + scenario = refScen.extractSingleScenario(scenNum); + assertTrue(isa(scenario, 'matRad_NominalScenario')); + ctScenIx = refScen.ctScenIx(scenNum); + ctScenNum = find(ctScenIx == refScen.ctScenProb(:, 1)); + assertEqual(scenario.ctScenProb, refScen.ctScenProb(ctScenNum, :)); + assertEqual(scenario.numOfCtScen, 1); + assertEqual(scenario.totNumScen, 1); + assertEqual(scenario.totNumShiftScen, 1); + assertEqual(scenario.totNumRangeScen, 1); + assertEqual(scenario.relRangeShift, refScen.relRangeShift(scenNum)); + assertEqual(scenario.absRangeShift, refScen.absRangeShift(scenNum)); + assertEqual(scenario.isoShift, refScen.isoShift(scenNum, :)); + assertEqual(scenario.maxAbsRangeShift, max(abs(refScen.absRangeShift(scenNum)))); + assertEqual(scenario.maxRelRangeShift, max(abs(refScen.relRangeShift(scenNum)))); + assertTrue(scenario.scenMask(ctScenIx, 1, 1)); + assertTrue(numel(find(scenario.scenMask)) == 1); + assertEqual(scenario.linearMask, [ctScenIx 1 1]); + assertElementsAlmostEqual(scenario.scenProb, helper_mvarGauss(scenario)); + assertEqual(scenario.scenForProb, refScen.scenForProb(scenNum, :)); + assertEqual(scenario.scenWeight, refScen.scenWeight(scenNum)); +end + function test_worstCaseScenarioCombineRange - model = matRad_WorstCaseScenarios(); +model = matRad_WorstCaseScenarios(); - assertExceptionThrown(@() helper_assignmentTest(model,'combineRange','hello'),'matRad:Error'); - assertTrue(model.combineRange); +assertExceptionThrown(@() helper_assignmentTest(model, 'combineRange', 'hello'), 'matRad:Error'); +assertTrue(model.combineRange); - assertEqual(model.totNumRangeScen,3); - model.combineRange = false; - assertFalse(model.combineRange); - assertEqual(model.totNumScen,15); - assertEqual(model.totNumRangeScen,9); +assertEqual(model.totNumRangeScen, 3); +model.combineRange = false; +assertFalse(model.combineRange); +assertEqual(model.totNumScen, 15); +assertEqual(model.totNumRangeScen, 9); function test_worstCaseScenarioShiftCombinations - model = matRad_WorstCaseScenarios(); - - assertExceptionThrown(@() helper_assignmentTest(model,'combinations','hello'),'matRad:Error'); - assertEqual(model.combinations,'none'); - - model.combinations = 'shift'; - assertEqual(model.combinations,'shift'); - assertEqual(model.totNumShiftScen,27); - assertEqual(model.totNumRangeScen,3); - assertEqual(model.totNumScen,29); - - model.combinations = 'all'; - assertEqual(model.combinations,'all'); - assertEqual(model.totNumShiftScen,27); - assertEqual(model.totNumRangeScen,3); - assertEqual(model.totNumScen,81); - - model.combineRange = false; - assertEqual(model.totNumShiftScen,27); - assertEqual(model.totNumRangeScen,9); - assertEqual(model.totNumScen,243); - - model.combinations = 'shift'; - assertEqual(model.totNumShiftScen,27); - assertEqual(model.totNumRangeScen,9); - assertEqual(model.totNumScen,35); \ No newline at end of file +model = matRad_WorstCaseScenarios(); + +assertExceptionThrown(@() helper_assignmentTest(model, 'combinations', 'hello'), 'matRad:Error'); +assertEqual(model.combinations, 'none'); + +model.combinations = 'shift'; +assertEqual(model.combinations, 'shift'); +assertEqual(model.totNumShiftScen, 27); +assertEqual(model.totNumRangeScen, 3); +assertEqual(model.totNumScen, 29); + +model.combinations = 'all'; +assertEqual(model.combinations, 'all'); +assertEqual(model.totNumShiftScen, 27); +assertEqual(model.totNumRangeScen, 3); +assertEqual(model.totNumScen, 81); + +model.combineRange = false; +assertEqual(model.totNumShiftScen, 27); +assertEqual(model.totNumRangeScen, 9); +assertEqual(model.totNumScen, 243); + +model.combinations = 'shift'; +assertEqual(model.totNumShiftScen, 27); +assertEqual(model.totNumRangeScen, 9); +assertEqual(model.totNumScen, 35); + +function test_worstCaseScenarioGridPoints +% Default: numOfSetupGridPoints=3, numOfRangeGridPoints=3 -> 9 scenarios +% (1 nominal + 6 shifts = 7 shift scenarios, 3 range scenarios) +model = matRad_WorstCaseScenarios(); +assertEqual(model.numOfSetupGridPoints, 3); +assertEqual(model.numOfRangeGridPoints, 3); +assertEqual(model.totNumShiftScen, 7); +assertEqual(model.totNumRangeScen, 3); +assertEqual(model.totNumScen, 9); + +% Invalid grid point values should throw an error +assertExceptionThrown(@() helper_assignmentTest(model, 'numOfSetupGridPoints', 2), 'matRad:Error'); +assertExceptionThrown(@() helper_assignmentTest(model, 'numOfRangeGridPoints', 2), 'matRad:Error'); + +% numOfRangeGridPoints=1: range worst cases removed -> 7 scenarios +% (7 shift scenarios, 1 range scenario) +model.numOfRangeGridPoints = 1; +assertEqual(model.numOfRangeGridPoints, 1); +assertEqual(model.totNumShiftScen, 7); +assertEqual(model.totNumRangeScen, 1); +assertEqual(model.totNumScen, 7); + +% numOfSetupGridPoints=1: shift worst cases removed -> 3 scenarios +% (1 shift scenario, 1 range scenario -> but range is still 1 from above) +% Reset range first, then set setup to 1 +model.numOfRangeGridPoints = 3; +model.numOfSetupGridPoints = 1; +assertEqual(model.numOfSetupGridPoints, 1); +assertEqual(model.totNumShiftScen, 1); +assertEqual(model.totNumRangeScen, 3); +assertEqual(model.totNumScen, 3); From 2e78813b8a8c7e12b1840ef7bf70089b87ab4610 Mon Sep 17 00:00:00 2001 From: Niklas Wahl Date: Tue, 31 Mar 2026 13:25:31 +0200 Subject: [PATCH 2/2] fix scenario index from subscript function and reformat files according to code style guidelines --- matRad/scenarios/matRad_ScenarioModel.m | 312 +++++++++++++----------- test/scenarios/test_nominalScenario.m | 158 ++++++------ test/scenarios/test_wcScenarios.m | 29 +++ 3 files changed, 279 insertions(+), 220 deletions(-) diff --git a/matRad/scenarios/matRad_ScenarioModel.m b/matRad/scenarios/matRad_ScenarioModel.m index c51f204d8..d96376a4a 100644 --- a/matRad/scenarios/matRad_ScenarioModel.m +++ b/matRad/scenarios/matRad_ScenarioModel.m @@ -1,75 +1,75 @@ classdef (Abstract) matRad_ScenarioModel < handle -% matRad_ScenarioModel -% This is an abstract interface class to define Scenario Models for use in -% robust treatment planning and uncertainty analysis. -% Subclasses should at least implement the update() function to generate -% their own scenarios. -% -% constructor (Abstract) -% matRad_ScenarioModel() -% matRad_ScenarioModel(ct) -% -% input: -% ct: ct cube -% -% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% -% Copyright 2022-2026 the matRad development team. -% -% This file is part of the matRad project. It is subject to the license -% terms in the LICENSE file found in the top-level directory of this -% distribution and at https://github.com/e0404/matRad/LICENSE.md. No part -% of the matRad project, including this file, may be copied, modified, -% propagated, or distributed except according to the terms contained in the -% LICENSE file. -% -% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - properties (AbortSet = true) %We use AbortSet = true here to avoid updates when - %Uncertainty model - rangeRelSD = 3.5; % given in % - rangeAbsSD = 1; % given in [mm] - shiftSD = [2.25 2.25 2.25]; % given in [mm] - wcSigma = 1; % Multiplier to compute the worst case / maximum shifts - - ctScenProb = [1 1]; % Ct Scenarios to be included in the model. Left column: Scenario Index. Right column: Scenario Probability + % matRad_ScenarioModel + % This is an abstract interface class to define Scenario Models for use in + % robust treatment planning and uncertainty analysis. + % Subclasses should at least implement the update() function to generate + % their own scenarios. + % + % constructor (Abstract) + % matRad_ScenarioModel() + % matRad_ScenarioModel(ct) + % + % input: + % ct: ct cube + % + % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % + % Copyright 2022-2026 the matRad development team. + % + % This file is part of the matRad project. It is subject to the license + % terms in the LICENSE file found in the top-level directory of this + % distribution and at https://github.com/e0404/matRad/LICENSE.md. No part + % of the matRad project, including this file, may be copied, modified, + % propagated, or distributed except according to the terms contained in the + % LICENSE file. + % + % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + properties (AbortSet = true) % We use AbortSet = true here to avoid updates when + % Uncertainty model + rangeRelSD = 3.5 % given in % + rangeAbsSD = 1 % given in [mm] + shiftSD = [2.25 2.25 2.25] % given in [mm] + wcSigma = 1 % Multiplier to compute the worst case / maximum shifts + + ctScenProb = [1 1] % Ct Scenarios to be included in the model. Left column: Scenario Index. Right column: Scenario Probability end - properties (Abstract,SetAccess = protected) + properties (Abstract, SetAccess = protected) name shortName end properties (Dependent) - wcFactor; + wcFactor end - + properties (SetAccess = protected) - numOfCtScen; % total number of CT scenarios used - numOfAvailableCtScen; % total number of CT scenarios existing in ct structure - ctScenIx; % map of all ct scenario indices per scenario - - - % these parameters will be filled according to the choosen scenario type - isoShift; - relRangeShift; - absRangeShift; - - maxAbsRangeShift; - maxRelRangeShift; - - totNumShiftScen; % total number of shift scenarios in x,y and z direction - totNumRangeScen; % total number of range and absolute range scenarios - totNumScen; % total number of samples - - scenForProb; % matrix for probability calculation - each row denotes one scenario, whereas columns denotes the realization value - scenProb; % probability of each scenario stored in a vector (according to uncertainty model) - scenWeight; % weight of scenario relative to the underlying uncertainty model (depends on how scenarios are chosen / sampled) - scenMask; - linearMask; + numOfCtScen % total number of CT scenarios used + numOfAvailableCtScen % total number of CT scenarios existing in ct structure + ctScenIx % map of all ct scenario indices per scenario + + % these parameters will be filled according to the chosen scenario type + isoShift + relRangeShift + absRangeShift + + maxAbsRangeShift + maxRelRangeShift + + totNumShiftScen % total number of shift scenarios in x,y and z direction + totNumRangeScen % total number of range and absolute range scenarios + totNumScen % total number of samples + + scenForProb % matrix for probability calculation - each row denotes one scenario, whereas columns denotes the realization value + scenProb % probability of each scenario stored in a vector (according to uncertainty model) + scenWeight % weight of scenario relative to the underlying uncertainty model (depends on how scenarios are chosen / sampled) + scenMask + linearMask end - + methods + function this = matRad_ScenarioModel(ct) if nargin == 0 || isempty(ct) this.numOfCtScen = 1; @@ -79,27 +79,28 @@ this.numOfAvailableCtScen = ct.numOfCtScen; end - this.ctScenProb = [(1:this.numOfCtScen)', ones(this.numOfCtScen,1)./this.numOfCtScen]; %Equal probability to be in each phase of the 4D ct - - %TODO: We could do this here automatically in the constructor, but - %Octave 5 has a bug here and throws an error - %this.updateScenarios(); + % Equal probability to be in each phase of the 4D ct + this.ctScenProb = [(1:this.numOfCtScen)', ones(this.numOfCtScen, 1) ./ this.numOfCtScen]; + + % TODO: We could do this here automatically in the constructor, but + % Octave 5 has a bug here and throws an error + % this.updateScenarios(); end function listAllScenarios(this) matRad_cfg = MatRad_Config.instance(); matRad_cfg.dispInfo('Listing all scenarios...\n'); matRad_cfg.dispInfo('\t#\tctScen\txShift\tyShift\tzShift\tabsRng\trelRng\tprob.\n'); - for s = 1:size(this.scenForProb,1) - str = [num2str(this.scenForProb(s,1),'%d\t'),sprintf('\t'), num2str(this.scenForProb(s,2:end),'\t%.3f')]; - matRad_cfg.dispInfo('\t%d\t%s\t%.3f\n',s,str,this.scenProb(s)); + for s = 1:size(this.scenForProb, 1) + str = [num2str(this.scenForProb(s, 1), '%d\t'), sprintf('\t'), num2str(this.scenForProb(s, 2:end), '\t%.3f')]; + matRad_cfg.dispInfo('\t%d\t%s\t%.3f\n', s, str, this.scenProb(s)); end end %% SETTERS & UPDATE - function set.rangeRelSD(this,rangeRelSD) + function set.rangeRelSD(this, rangeRelSD) valid = isnumeric(rangeRelSD) && isscalar(rangeRelSD) && rangeRelSD >= 0; - if ~valid + if ~valid matRad_cfg = MatRad_Config.instance(); matRad_cfg.dispError('Invalid value for rangeRelSD! Needs to be a real positive scalar!'); end @@ -107,9 +108,9 @@ function listAllScenarios(this) this.updateScenarios(); end - function set.rangeAbsSD(this,rangeAbsSD) + function set.rangeAbsSD(this, rangeAbsSD) valid = isnumeric(rangeAbsSD) && isscalar(rangeAbsSD) && rangeAbsSD >= 0; - if ~valid + if ~valid matRad_cfg = MatRad_Config.instance(); matRad_cfg.dispError('Invalid value for rangeAbsSD! Needs to be a real positive scalar!'); end @@ -117,9 +118,9 @@ function listAllScenarios(this) this.updateScenarios(); end - function set.shiftSD(this,shiftSD) + function set.shiftSD(this, shiftSD) valid = isnumeric(shiftSD) && isrow(shiftSD) && numel(shiftSD) == 3 && all(shiftSD > 0); - if ~valid + if ~valid matRad_cfg = MatRad_Config.instance(); matRad_cfg.dispError('Invalid value for shiftSD! Needs to be 3-element numeric row vector!'); end @@ -127,9 +128,9 @@ function listAllScenarios(this) this.updateScenarios(); end - function set.wcSigma(this,wcSigma) + function set.wcSigma(this, wcSigma) valid = isnumeric(wcSigma) && isscalar(wcSigma) && wcSigma >= 0; - if ~valid + if ~valid matRad_cfg = MatRad_Config.instance(); matRad_cfg.dispError('Invalid value for wcSigma! Needs to be a real positive scalar!'); end @@ -137,70 +138,83 @@ function listAllScenarios(this) this.updateScenarios(); end - function set.ctScenProb(this,ctScenProb) - valid = isnumeric(ctScenProb) && ismatrix(ctScenProb) && size(ctScenProb,2) == 2 && all(round(ctScenProb(:,1)) == ctScenProb(:,1)) && all(ctScenProb(:) >= 0); + function set.ctScenProb(this, ctScenProb) + valid = isnumeric(ctScenProb) && ... + ismatrix(ctScenProb) && ... + size(ctScenProb, 2) == 2 && ... + all(round(ctScenProb(:, 1)) == ctScenProb(:, 1)) && ... + all(ctScenProb(:) >= 0); + if ~valid matRad_cfg = MatRad_Config.instance(); - matRad_cfg.dispError('Invalid value for used ctScenProb! Needs to be a valid 2-column matrix with left column representing the scenario index and right column representing the appropriate probabilities [0,1]!'); - end + matRad_cfg.dispError(['Invalid value for used ctScenProb! ' ... + 'Needs to be a valid 2-column matrix with left column representing the scenario index ' ... + 'and right column representing the appropriate probabilities [0,1]!']); + end this.ctScenProb = ctScenProb; this.updateScenarios(); end - - function scenarios = updateScenarios(this) - %This function will always update the scenarios given the - %current property settings + function scenarios = updateScenarios(this) + % This function will always update the scenarios given the + % current property settings matRad_cfg = MatRad_Config.instance(); matRad_cfg.dispError('This abstract function needs to be implemented!'); end - function newInstance = extractSingleScenario(this,scenNum) + function newInstance = extractSingleScenario(this, scenNum) newInstance = matRad_NominalScenario(); - - ctScenNum = this.linearMask(scenNum,1); - - %First set properties that force an update - newInstance.numOfCtScen = 1; - newInstance.ctScenProb = this.ctScenProb(ctScenNum,:); - - %Now overwrite existing variables for correct probabilties and - %error realizations - newInstance.scenForProb = this.scenForProb(scenNum,:); - newInstance.relRangeShift = this.scenForProb(scenNum,6); - newInstance.absRangeShift = this.scenForProb(scenNum,5); - newInstance.isoShift = this.scenForProb(scenNum,2:4); + + ctScenNum = this.linearMask(scenNum, 1); + + % First set properties that force an update + newInstance.numOfCtScen = 1; + newInstance.ctScenProb = this.ctScenProb(ctScenNum, :); + + % Now overwrite existing variables for correct probabilities + % and error realizations + newInstance.scenForProb = this.scenForProb(scenNum, :); + newInstance.relRangeShift = this.scenForProb(scenNum, 6); + newInstance.absRangeShift = this.scenForProb(scenNum, 5); + newInstance.isoShift = this.scenForProb(scenNum, 2:4); newInstance.scenProb = this.scenProb(scenNum); newInstance.scenWeight = this.scenWeight(scenNum); newInstance.maxAbsRangeShift = max(abs(this.absRangeShift(scenNum))); newInstance.maxRelRangeShift = max(abs(this.relRangeShift(scenNum))); - newInstance.scenMask = false(this.numOfAvailableCtScen,1,1); + newInstance.scenMask = false(this.numOfAvailableCtScen, 1, 1); newInstance.linearMask = [newInstance.ctScenIx 1 1]; - - newInstance.scenMask(newInstance.linearMask(:,1),newInstance.linearMask(:,2),newInstance.linearMask(:,3)) = true; - %newInstance.updateScenarios(); + + newInstance.scenMask(newInstance.linearMask(:, 1), newInstance.linearMask(:, 2), newInstance.linearMask(:, 3)) = true; + % newInstance.updateScenarios(); end - - function scenIx = sub2scenIx(this,ctScen,shiftScen,rangeShiftScen) - %Returns linear index in the scenario cell array from scenario - %subscript indices - if ~isvector(this.scenMask) - scenIx = sub2ind(size(this.scenMask),ctScen,shiftScen,rangeShiftScen); - else + + function scenIx = sub2scenIx(this, ctScen, shiftScen, rangeShiftScen) + % Returns linear index in the scenario cell array from scenario + % subscript indices. + % Note: + % scenMask can be squeezed by MATLAB (e.g. [1 x N_S x 1] + % becomes [1 x N_S]), so we must explicitly check for a + % column vector here in case of 4D scenarios only + if iscolumn(this.scenMask) + % pure CT scenarios scenIx = this.ctScenIx(ctScen); + else + % multi-dimensional scenario array + scenIx = sub2ind(size(this.scenMask), ctScen, shiftScen, rangeShiftScen); end end - function scenNum = scenNum(this,fullScenIx) - %gets number of scneario from full scenario index in scenMask + function scenNum = scenNum(this, fullScenIx) + % gets number of scenario from full scenario index in scenMask scenNum = find(find(this.scenMask) == fullScenIx); end - + %% Deprecated functions / properties - function newInstance = extractSingleNomScen(this,~,scenIdx) + function newInstance = extractSingleNomScen(this, ~, scenIdx) matRad_cfg = MatRad_Config.instance(); - matRad_cfg.dispDeprecationWarning('The function extractSingleNomScen of the scenario class will soon be deprecated! Use extractSingleScenario instead!'); + matRad_cfg.dispDeprecationWarning(['The function extractSingleNomScen of the scenario class will soon be deprecated! ' ... + 'Use extractSingleScenario instead!']); newInstance = this.extractSingleScenario(scenIdx); end @@ -216,7 +230,7 @@ function listAllScenarios(this) value = this.wcSigma; end - function set.wcFactor(this,value) + function set.wcFactor(this, value) matRad_cfg = MatRad_Config.instance(); matRad_cfg.dispDeprecationWarning('The property wcFactor of the scenario class will soon be deprecated!'); this.wcSigma = value; @@ -225,59 +239,59 @@ function listAllScenarios(this) end methods (Static) - + function classList = getAvailableModels() matRad_cfg = MatRad_Config.instance(); - - %Use the root folder and the scenarios folder only + + % Use the root folder and the scenarios folder only folders = {fileparts(mfilename('fullpath'))}; folders = [folders matRad_cfg.userfolders]; persistent metaScenarioModels lastOptionalPaths - - %First we do a sanity check if persistently stored metaclasses are valid - if ~matRad_cfg.isOctave && ~isempty(metaScenarioModels) && ~all(cellfun(@isvalid,metaScenarioModels)) + + % First we do a sanity check if persistently stored metaclasses are valid + if ~matRad_cfg.isOctave && ~isempty(metaScenarioModels) && ~all(cellfun(@isvalid, metaScenarioModels)) matRad_cfg.dispWarning('Found invalid ScenarioModels, updating model cache.'); metaScenarioModels = []; end if isempty(metaScenarioModels) || (~isempty(lastOptionalPaths) && ~isequal(lastOptionalPaths, folders)) lastOptionalPaths = folders; - metaScenarioModels = matRad_findSubclasses(meta.class.fromName(mfilename('class')),'folders',folders,'includeSubfolders',true); + metaScenarioModels = matRad_findSubclasses(meta.class.fromName(mfilename('class')), 'folders', folders, 'includeSubfolders', true); end - classList = matRad_identifyClassesByConstantProperties(metaScenarioModels,'shortName','defaults',{'nomScen'}); + classList = matRad_identifyClassesByConstantProperties(metaScenarioModels, 'shortName', 'defaults', {'nomScen'}); if isempty(classList) - matRad_cfg.dispError('No models found in paths %s',strjoin(folders,'\n')); + matRad_cfg.dispError('No models found in paths %s', strjoin(folders, '\n')); end end - - function model = create(modelMetadata,ct) - if isa(modelMetadata,'matRad_ScenarioModel') + + function model = create(modelMetadata, ct) + if isa(modelMetadata, 'matRad_ScenarioModel') model = modelMetadata; - return; + return end - + matRad_cfg = MatRad_Config.instance(); if ischar(modelMetadata) || isstring(modelMetadata) - modelMetadata = struct('model',modelMetadata); + modelMetadata = struct('model', modelMetadata); end modelClassList = matRad_ScenarioModel.getAvailableModels(); modelNames = {modelClassList.shortName}; - - if ~isfield(modelMetadata,'model') || ~any(strcmp(modelNames,modelMetadata.model)) + + if ~isfield(modelMetadata, 'model') || ~any(strcmp(modelNames, modelMetadata.model)) matRad_cfg.dispWarning('Scenario Model not found, creating nominal scenario instead!'); modelMetadata.model = 'nomScen'; end - usedModel = find(strcmp(modelNames,modelMetadata.model)); - + usedModel = find(strcmp(modelNames, modelMetadata.model)); + if ~isscalar(usedModel) usedModel = usedModel(1); end - + modelClassInfo = modelClassList(usedModel); if nargin < 2 @@ -286,24 +300,24 @@ function listAllScenarios(this) model = modelClassInfo.handle(ct); end - modelMetadata = rmfield(modelMetadata,'model'); - - %Now overwrite properties + modelMetadata = rmfield(modelMetadata, 'model'); + + % Now overwrite properties fields = fieldnames(modelMetadata); - + % iterate over all fieldnames and try to set the % corresponding properties inside the engine if matRad_cfg.isOctave - c2sWarningState = warning('off','Octave:classdef-to-struct'); + c2sWarningState = warning('off', 'Octave:classdef-to-struct'); end - + for i = 1:length(fields) try field = fields{i}; - if matRad_ispropCompat(model,field) - model.(field) = matRad_recursiveFieldAssignment(model.(field),modelMetadata.(field),true); + if matRad_ispropCompat(model, field) + model.(field) = matRad_recursiveFieldAssignment(model.(field), modelMetadata.(field), true); else - matRad_cfg.dispWarning('Not able to assign property ''%s'' from multScen struct to Scenario Model!',field); + matRad_cfg.dispWarning('Not able to assign property ''%s'' from multScen struct to Scenario Model!', field); end catch ME % catch exceptions when the model has no properties, @@ -314,17 +328,17 @@ function listAllScenarios(this) matRad_cfg = MatRad_Config.instance(); switch ME.identifier case 'MATLAB:noPublicFieldForClass' - matRad_cfg.dispWarning('Not able to assign property from multScen struct to scenario model: %s',ME.message); + matRad_cfg.dispWarning('Not able to assign property from multScen struct to scenario model: %s', ME.message); otherwise - matRad_cfg.dispWarning('Problem while setting up scenario Model from struct:%s %s',field,ME.message); + matRad_cfg.dispWarning('Problem while setting up scenario Model from struct:%s %s', field, ME.message); end end end - + if matRad_cfg.isOctave - warning(c2sWarningState.state,'Octave:classdef-to-struct'); + warning(c2sWarningState.state, 'Octave:classdef-to-struct'); end end + end end - diff --git a/test/scenarios/test_nominalScenario.m b/test/scenarios/test_nominalScenario.m index 10e866a31..60d685af2 100644 --- a/test/scenarios/test_nominalScenario.m +++ b/test/scenarios/test_nominalScenario.m @@ -1,84 +1,100 @@ function test_suite = test_nominalScenario -test_functions=localfunctions(); +test_functions = localfunctions(); initTestSuite; function test_nominalScenarioConstructor - scenario = matRad_NominalScenario(); - assertTrue(isa(scenario, 'matRad_NominalScenario')); - assertTrue(isa(scenario, 'matRad_ScenarioModel')); - assertEqual(scenario.shortName, 'nomScen'); - assertEqual(scenario.ctScenProb, [1 1]); - assertEqual(scenario.numOfCtScen, 1); - assertEqual(scenario.totNumScen, 1); - assertEqual(scenario.totNumShiftScen, 1); - assertEqual(scenario.totNumRangeScen, 1); - assertEqual(scenario.relRangeShift, 0); - assertEqual(scenario.absRangeShift, 0); - assertEqual(scenario.isoShift, zeros(1,3)); - assertEqual(scenario.maxAbsRangeShift, 0); - assertEqual(scenario.maxRelRangeShift, 0); - assertEqual(scenario.scenMask, true(1,1,1)); - assertEqual(scenario.linearMask, [1 1 1]); - assertElementsAlmostEqual(scenario.scenProb,helper_mvarGauss(scenario)); - assertEqual(scenario.scenForProb,[1 zeros(1,5)]); - assertEqual(scenario.scenWeight, 1); +scenario = matRad_NominalScenario(); +assertTrue(isa(scenario, 'matRad_NominalScenario')); +assertTrue(isa(scenario, 'matRad_ScenarioModel')); +assertEqual(scenario.shortName, 'nomScen'); +assertEqual(scenario.ctScenProb, [1 1]); +assertEqual(scenario.numOfCtScen, 1); +assertEqual(scenario.totNumScen, 1); +assertEqual(scenario.totNumShiftScen, 1); +assertEqual(scenario.totNumRangeScen, 1); +assertEqual(scenario.relRangeShift, 0); +assertEqual(scenario.absRangeShift, 0); +assertEqual(scenario.isoShift, zeros(1, 3)); +assertEqual(scenario.maxAbsRangeShift, 0); +assertEqual(scenario.maxRelRangeShift, 0); +assertEqual(scenario.scenMask, true(1, 1, 1)); +assertEqual(scenario.linearMask, [1 1 1]); +assertElementsAlmostEqual(scenario.scenProb, helper_mvarGauss(scenario)); +assertEqual(scenario.scenForProb, [1 zeros(1, 5)]); +assertEqual(scenario.scenWeight, 1); function test_nominalScenarioConstructorWithCt - n = 5; - ct = struct('numOfCtScen',n); - scenario = matRad_NominalScenario(ct); - assertTrue(isa(scenario, 'matRad_NominalScenario')); - assertTrue(isa(scenario, 'matRad_ScenarioModel')); - assertEqual(scenario.shortName, 'nomScen'); - assertEqual(scenario.ctScenProb, [(1:n)' ones(n,1)./n]); - assertEqual(scenario.numOfCtScen, n); - assertEqual(scenario.totNumScen, n); - assertEqual(scenario.totNumShiftScen, 1); - assertEqual(scenario.totNumRangeScen, 1); - assertEqual(scenario.relRangeShift, zeros(n,1)); - assertEqual(scenario.absRangeShift, zeros(n,1)); - assertEqual(scenario.isoShift, zeros(n,3)); - assertEqual(scenario.maxAbsRangeShift, 0); - assertEqual(scenario.maxRelRangeShift, 0); - assertTrue(isequal(scenario.scenMask, true(n,1,1))); - assertEqual(scenario.linearMask, [(1:n)' ones(n,2)]); - assertElementsAlmostEqual(scenario.scenProb,helper_mvarGauss(scenario)); - assertEqual(scenario.scenForProb,[(1:n)' zeros(n,5)]); - assertEqual(scenario.scenWeight, ones(n,1)./n); +n = 5; +ct = struct('numOfCtScen', n); +scenario = matRad_NominalScenario(ct); +assertTrue(isa(scenario, 'matRad_NominalScenario')); +assertTrue(isa(scenario, 'matRad_ScenarioModel')); +assertEqual(scenario.shortName, 'nomScen'); +assertEqual(scenario.ctScenProb, [(1:n)' ones(n, 1) ./ n]); +assertEqual(scenario.numOfCtScen, n); +assertEqual(scenario.totNumScen, n); +assertEqual(scenario.totNumShiftScen, 1); +assertEqual(scenario.totNumRangeScen, 1); +assertEqual(scenario.relRangeShift, zeros(n, 1)); +assertEqual(scenario.absRangeShift, zeros(n, 1)); +assertEqual(scenario.isoShift, zeros(n, 3)); +assertEqual(scenario.maxAbsRangeShift, 0); +assertEqual(scenario.maxRelRangeShift, 0); +assertTrue(isequal(scenario.scenMask, true(n, 1, 1))); +assertEqual(scenario.linearMask, [(1:n)' ones(n, 2)]); +assertElementsAlmostEqual(scenario.scenProb, helper_mvarGauss(scenario)); +assertEqual(scenario.scenForProb, [(1:n)' zeros(n, 5)]); +assertEqual(scenario.scenWeight, ones(n, 1) ./ n); -function test_nominalScenarioExtractSingleScenario - scenario = matRad_NominalScenario(); - newInstance = scenario.extractSingleScenario(1); - assertTrue(isequal(scenario,newInstance)); +function test_nominalScenarioSub2scenIx +% Nominal single scenario: scenMask is scalar true(1,1,1). +% sub2scenIx must return 1 without error. +scenario = matRad_NominalScenario(); +scenIx = scenario.sub2scenIx(1, 1, 1); +assertEqual(scenIx, 1); +assertEqual(scenario.scenNum(scenIx), 1); -function test_nominalScenarioExtractSingleScenarioWithCtScen - n = 5; - ct = struct('numOfCtScen',n); - refScen = matRad_NominalScenario(ct); +% Nominal with multiple CT: scenMask is true(n,1,1), squeezed to [n x 1] +% (column vector). iscolumn branch must route to ctScenIx correctly. +n = 5; +ct = struct('numOfCtScen', n); +scenarioCt = matRad_NominalScenario(ct); +for s = 1:scenarioCt.totNumScen + scenIx = scenarioCt.sub2scenIx(scenarioCt.linearMask(s, 1), scenarioCt.linearMask(s, 2), scenarioCt.linearMask(s, 3)); + assertEqual(scenarioCt.scenNum(scenIx), s); +end - for scenNum = 1:refScen.totNumScen - scenario = refScen.extractSingleScenario(scenNum); - assertTrue(isa(scenario, 'matRad_NominalScenario')); - ctScenIx = refScen.ctScenIx(scenNum); - ctScenNum = find(ctScenIx == refScen.ctScenProb(:,1)); - assertEqual(scenario.ctScenProb, refScen.ctScenProb(ctScenNum,:)); - assertEqual(scenario.numOfCtScen, 1); - assertEqual(scenario.totNumScen, 1); - assertEqual(scenario.totNumShiftScen, 1); - assertEqual(scenario.totNumRangeScen, 1); - assertEqual(scenario.relRangeShift, refScen.relRangeShift(scenNum)); - assertEqual(scenario.absRangeShift, refScen.absRangeShift(scenNum)); - assertEqual(scenario.isoShift, refScen.isoShift(scenNum,:)); - assertEqual(scenario.maxAbsRangeShift, max(abs(refScen.absRangeShift(scenNum)))); - assertEqual(scenario.maxRelRangeShift, max(abs(refScen.relRangeShift(scenNum)))); - assertTrue(scenario.scenMask(ctScenIx,1,1)); - assertTrue(numel(find(scenario.scenMask)) == 1); - assertEqual(scenario.linearMask, [ctScenIx 1 1]); - assertElementsAlmostEqual(scenario.scenProb,helper_mvarGauss(scenario)); - assertEqual(scenario.scenForProb,refScen.scenForProb(scenNum,:)); - assertEqual(scenario.scenWeight, refScen.scenWeight(scenNum)); - end +function test_nominalScenarioExtractSingleScenario +scenario = matRad_NominalScenario(); +newInstance = scenario.extractSingleScenario(1); +assertTrue(isequal(scenario, newInstance)); +function test_nominalScenarioExtractSingleScenarioWithCtScen +n = 5; +ct = struct('numOfCtScen', n); +refScen = matRad_NominalScenario(ct); +for scenNum = 1:refScen.totNumScen + scenario = refScen.extractSingleScenario(scenNum); + assertTrue(isa(scenario, 'matRad_NominalScenario')); + ctScenIx = refScen.ctScenIx(scenNum); + ctScenNum = find(ctScenIx == refScen.ctScenProb(:, 1)); + assertEqual(scenario.ctScenProb, refScen.ctScenProb(ctScenNum, :)); + assertEqual(scenario.numOfCtScen, 1); + assertEqual(scenario.totNumScen, 1); + assertEqual(scenario.totNumShiftScen, 1); + assertEqual(scenario.totNumRangeScen, 1); + assertEqual(scenario.relRangeShift, refScen.relRangeShift(scenNum)); + assertEqual(scenario.absRangeShift, refScen.absRangeShift(scenNum)); + assertEqual(scenario.isoShift, refScen.isoShift(scenNum, :)); + assertEqual(scenario.maxAbsRangeShift, max(abs(refScen.absRangeShift(scenNum)))); + assertEqual(scenario.maxRelRangeShift, max(abs(refScen.relRangeShift(scenNum)))); + assertTrue(scenario.scenMask(ctScenIx, 1, 1)); + assertTrue(numel(find(scenario.scenMask)) == 1); + assertEqual(scenario.linearMask, [ctScenIx 1 1]); + assertElementsAlmostEqual(scenario.scenProb, helper_mvarGauss(scenario)); + assertEqual(scenario.scenForProb, refScen.scenForProb(scenNum, :)); + assertEqual(scenario.scenWeight, refScen.scenWeight(scenNum)); +end diff --git a/test/scenarios/test_wcScenarios.m b/test/scenarios/test_wcScenarios.m index 3d646db99..67ae11b23 100644 --- a/test/scenarios/test_wcScenarios.m +++ b/test/scenarios/test_wcScenarios.m @@ -206,3 +206,32 @@ assertEqual(model.totNumShiftScen, 1); assertEqual(model.totNumRangeScen, 3); assertEqual(model.totNumScen, 3); + +function test_worstCaseScenarioSub2scenIxDegenerate +% Regression test for sub2scenIx when scenMask dimensions are degenerate +% (i.e., MATLAB squeezes trailing/leading singletons). + +% Case 1: numOfRangeGridPoints=1 -> scenMask is [1 x N_S x 1], MATLAB +% squeezes to [1 x N_S] (row vector). sub2scenIx must NOT fall into the +% CT-only branch (iscolumn = false for row vectors). +model = matRad_WorstCaseScenarios(); +model.numOfRangeGridPoints = 1; +assertEqual(model.totNumRangeScen, 1); +assertEqual(model.totNumShiftScen, 7); +assertEqual(model.numOfCtScen, 1); +% Verify sub2scenIx round-trips correctly for every scenario +for s = 1:model.totNumScen + scenIx = model.sub2scenIx(model.linearMask(s, 1), model.linearMask(s, 2), model.linearMask(s, 3)); + assertEqual(model.scenNum(scenIx), s); +end + +% Case 2: numOfSetupGridPoints=1 -> scenMask is [1 x 1 x N_R], stays 3-D. +% sub2ind path must work with the known dimensions, not size(scenMask). +model2 = matRad_WorstCaseScenarios(); +model2.numOfSetupGridPoints = 1; +assertEqual(model2.totNumShiftScen, 1); +assertEqual(model2.totNumRangeScen, 3); +for s = 1:model2.totNumScen + scenIx = model2.sub2scenIx(model2.linearMask(s, 1), model2.linearMask(s, 2), model2.linearMask(s, 3)); + assertEqual(model2.scenNum(scenIx), s); +end