Started implementing robust methodology. Framework completed.

NLP_beamlet_optimizer.m

@@ -1,13 +1,11 @@
 
-% this is a simple example for a constrained NLP solver from matlabs
-% optimization toolbox
-
 function [D_full, w_fin] = NLP_beamlet_optimizer
 % This function performs the beamlet optimization
 % Inputs: none. Still needs to find "Geometry" and "beamlets" from Wiscplan
 % Outputs: full dose image dose: D_full, optimal beamlet weights: w_fin
 %
-% Made by Peter Ferjancic May 2nd 2018
+% Made by Peter Ferjancic 1. May 2018
+% Last updated: 1. May 2018
 % Inspired by Ana Barrigan's REGGUI optimization procedures
 
 % To-do:
@@ -52,8 +50,13 @@ optGoal = make_ROI_goals(Geometry, beamlets);
 w0 = ones(numBeam,1);
 w0 = mean(optGoal{1}.target ./ (optGoal{1}.beamlets_pruned*w0)) * w0;
 
+% -- MAKE IT ROBUST --
+RO_params=0;
+optGoal = make_robust_optGoal(optGoal, RO_params);
+
+
 % -- CALLBACK OPTIMIZATION FUNCTION --
-fun = @(x) eval_f(x, optGoal);
+fun = @(x) get_penalty(x, optGoal);
 
 % -- OPTIMIZATION PARAMETERS --
 % define optimization parameters
@@ -67,7 +70,7 @@ nonlcon = [];
 
 % define opt limits, and make it fmincon progress
 options = optimoptions('fmincon');
-options.MaxFunctionEvaluations = 100000;
+options.MaxFunctionEvaluations = 30000;
 options.Display = 'iter';
 options.PlotFcn = 'optimplotfval';
 fprintf('-')
@@ -95,27 +98,55 @@ end
 % orthoslice(D_full, [0,70])
 
 %% support functions
-% ---- PENALTY FUNCTION ----
-function penalty = eval_f(x, optGoal)
+% ---- PENALTY FUNCTION 1 ----
+function penalty = get_penalty(x, optGoal)
+    % this function gets called by the optimizer. It checks the penalty for
+    % all the robust implementation and returns the worst result.
+    
+    NumScenarios = optGoal{1}.NbrRandScenarios * optGoal{1}.NbrSystSetUpScenarios * optGoal{1}.NbrRangeScenarios;
+    fobj = zeros(NumScenarios,1);  
+    sc_i = 1;
+    
+    for nrs_i = 1:optGoal{1}.NbrRandScenarios 
+        for sss_i = 1:optGoal{1}.NbrSystSetUpScenarios % syst. setup scenarios = ss
+            for rrs_i = 1:optGoal{1}.NbrRangeScenarios % range scenario = rs
+                fobj(sc_i)=eval_f(x, optGoal, nrs_i, sss_i, rrs_i);
+                sc_i = sc_i + 1;
+            end
+        end
+    end
+    % take the worst case penalty of evaluated scenarios
+    penalty=max(fobj);
+end
+
+% ---- Evaluate penalty for single scenario ----
+function penalty = eval_f(x, optGoal, nrs_i, sss_i, rrs_i)
     penalty = 0;
     % for each condition
-    for i = 1:numel(optGoal)
-        switch optGoal{i}.function
+    for goal_i = 1:numel(optGoal)
+        switch optGoal{goal_i}.function
             case 'min'
                 % penalize if achieved dose is lower than target dose
-                d_penalty = sum(max(0, optGoal{i}.target-(optGoal{i}.beamlets_pruned * x)));
+                d_penalty = 1.0e0 * sum(max(0, ...
+                    (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.target) -...
+                    (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.beamlets_pruned * x)));
             case 'max'
                 % penalize if achieved dose is higher than target dose
-                d_penalty = sum(max(0, (optGoal{i}.beamlets_pruned * x)-(optGoal{i}.target+5)));
+                d_penalty = 1.0e0 * sum(max(0, ...
+                    (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.beamlets_pruned * x)-...
+                    (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.target)));
             case 'LeastSquare'
                 % penalize with sum of squares any deviation from target
                 % dose
-                d_penalty = sum(((beamlets * x) - target).^2);
+                d_penalty = 1.0e-1* sum(((...
+                    optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.beamlets_pruned * x) - ...
+                    optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.target).^2);
         end
-        penalty = penalty + d_penalty * optGoal{i}.opt_weight;
+        penalty = penalty + d_penalty * optGoal{goal_i}.opt_weight;
     end
 end
 
+
 % ---- DVH PLOT FUNCTION ----
 function plot_DVH(dose, optGoal, optGoal_idx)
     % this function plots the DVHs of the given dose
@@ -188,7 +219,7 @@ function optGoal = make_ROI_goals(Geometry, beamlets)
     ROI_idx = Geometry.ROIS{1, 3}.ind;
     goal_3.ROI_idx = ROI_idx;
     goal_3.D_final = 0;
-    goal_3.function = 'max';
+    goal_3.function = 'LeastSquare';
     goal_3.beamlets_pruned = beamlets(ROI_idx, :);
     goal_3.target   = ones(numel(ROI_idx), 1) * goal_3.D_final;
     goal_3.opt_weight = 1 / numel(ROI_idx); % normalize to volume of target area
@@ -198,6 +229,63 @@ function optGoal = make_ROI_goals(Geometry, beamlets)
     % -- END DEFINITION OF GOAL --
     
 end
+% ---- MAKE ROI ROBUST ----
+function optGoal = make_robust_optGoal(optGoal, RO_params);
+    % take regular optimal goal and translate it into several robust cases
+    
+    % nrs - random scenarios
+    % sss - system setup scenarios
+    % rrs - random range scenarios
+    nrs_scene_list={[0,0,0]};
+    sss_scene_list={[0,0,0], [5,5,5], [-5, -5, -5]};
+    rrs_scene_list={[0,0,0]};
+    
+    for i = 1:numel(optGoal)
+        optGoal{i}.NbrRandScenarios     =numel(nrs_scene_list);
+        optGoal{i}.NbrSystSetUpScenarios=numel(sss_scene_list);
+        optGoal{i}.NbrRangeScenarios    =numel(rrs_scene_list);
+    end
+    
 
+    for goal_i = 1:numel(optGoal)
+        out_01.beamlets_pruned  = optGoal{goal_i}.beamlets_pruned;
+        out_01.target           = optGoal{goal_i}.target;
+        
+        for nrs_i = 1:optGoal{goal_i}.NbrRandScenarios          % num. of random scenarios
+            
+            out_02.beamlets_pruned  = out_01.beamlets_pruned;
+            out_02.target           = out_01.target;
+            
+            for sss_i = 1:optGoal{goal_i}.NbrSystSetUpScenarios   % syst. setup scenarios = sss
+                
+                out_03 = get_RO_sss(out_02);
+                
+                for rrs_i = 1:optGoal{goal_i}.NbrRangeScenarios   % range scenario = rrs
+                    
+                    out_final.beamlets_pruned = out_03.beamlets_pruned;
+                    out_final.target = out_03.target;
+                    
+                    optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.beamlets_pruned   = optGoal{goal_i}.beamlets_pruned;
+                    optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.target     = optGoal{goal_i}.target;
+                end
+            end
+        end
+    end
+end
 
+function out_03 = get_RO_sss(out_02);
+    % this function returns a systemic perturbation of the original data
+    
+    %% THIS PART DOES NOT WORK WELL! NEEDS TO BE FIXED SO THAT IT ACTUALLY PERMUTES DATA
+%     wbar1 = waitbar(0, 'Creating SSS beamlet array');
+%     for beam_i=1:numBeam
+%         % for each beam define how much dose it delivers on each voxel
+%         idx=beamlet_cell_array{1, beam_i}.non_zero_indices;
+%         beamlets(idx, beam_i) = 1000*beamlet_cell_array{1, beam_i}.non_zero_values;
+%         waitbar(beam_i/numBeam)
+%     end
+%     close(wbar1)
+
+    out_03= out_02;
+end