WiscPlan_v2/WiscPlanPhotonkV125/matlab_frontend/NLP_beamlet_optimizer.m

% colorwash(Geometry.data, D_full, [500, 1500], [0,70])
% orthoslice(D_full, [0,70])

function [D_full, w_fin, Geometry, optGoal] = 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
% 
% [D_full, w_fin, Geometry, optGoal] = NLP_beamlet_optimizer;
% 
% Made by Peter Ferjancic 1. May 2018
% Last updated: 14. August 2018
% Inspired by Ana Barrigan's REGGUI optimization procedures

N_fcallback1 = 5000;
N_fcallback2 = 200000;
patient = 'temp';
goalsName= 'ROI_goals_DP';

switch patient
    case 'temp'
        patient_dir = 'C:\000-temp\DELETEME';
        DP_dir = '\\Mpufs5\data_wnx1\_Data\Avastin\AV009\PF_RODP_analysis'
    case 'avastin_009'
        patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatData_Avastin009';
        DP_dir = '\\Mpufs5\data_wnx1\_Data\Avastin\AV009\PF_RODP_analysis';
    case 'gbm_005'
%         patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatData_ausGli_005';
        patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatData_ausGli_005_16beam_2';
        DP_dir = '\\Mpufs5\data_wnx1\_Data\Glioma_aus\FET_FGL005\B1\Processed';
    case 'gbm_015'
        patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatData_ausGli_015';
        DP_dir = '\\Mpufs5\data_wnx1\_Data\Glioma_aus\FET_FGL015\B1\Processed';
    case 'gbm_022'
        patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatData_ausGli_022';
        DP_dir = '\\Mpufs5\data_wnx1\_Data\Glioma_aus\FET_FGL022\B1\Processed';
    case 'avastin_009'
    otherwise
        error('invalid case')
end

% merge_beamlets(4, patient_dir);


%% PROGRAM STARTS HERE
% - no tocar lo que hay debajo -
fprintf('starting NLP optimization process... \n')

% % -- LOAD GEOMETRY AND BEAMLETS --
load([patient_dir '\matlab_files\Geometry.mat'])


%% -- OPTIMIZATION TARGETS --
load([DP_dir '\RODP_files\' goalsName '.mat']);

optGoal = ROI_goals.optGoal;
optGoal_beam = ROI_goals.optGoal_beam;
optGoal_idx = ROI_goals.optGoal_idx;
targetMinMax_idx = ROI_goals.targetMinMax_idx;

[beamlets, beamlets_joined, numBeamlet, numBeam, beam_i_list] = get_beam_lets(Geometry, patient_dir);


% -- make them robust --
RO_params=0;
optGoal_beam = make_robust_optGoal(optGoal_beam, RO_params, beamlets_joined);
optGoal = make_robust_optGoal(optGoal, RO_params, beamlets);

%% -- INITIALIZE BEAMLET WEIGHTS --
w0_beams = ones(numBeam,1);
w0_beams = mean(optGoal_beam{1}.target ./ (optGoal_beam{1}.beamlets_pruned*w0_beams+0.1)) * w0_beams;
w0_beams = w0_beams + (2*rand(size(w0_beams))-1) *0.1 .*w0_beams; % random perturbation
w0_beams = double(w0_beams);

% -- CALLBACK OPTIMIZATION FUNCTION --
fun1 = @(x) get_penalty(x, optGoal_beam);
fun2 = @(x) get_penalty(x, optGoal);

% -- OPTIMIZATION PARAMETERS --
% define optimization parameters
A = [];
b = [];
Aeq = [];
beq = [];
lb = zeros(1, numBeamlet);
lb_beam = zeros(1, numBeam);
ub = [];
nonlcon = [];

% define opt limits, and make it fmincon progress
options = optimoptions('fmincon');
options.MaxFunctionEvaluations = N_fcallback1;
options.Display = 'iter';
options.PlotFcn = 'optimplotfval';
% options.UseParallel = true;
% options.OptimalityTolerance = 1e-9;

fprintf('\n running initial optimizer:')

%% Run the optimization
% -- GET FULL BEAM WEIGHTS --
tic
w_beam = fmincon(fun1,w0_beams,A,b,Aeq,beq,lb_beam,ub,nonlcon,options);
fprintf('  done!:')
% t=toc;
% disp(['Optimization time for beams = ',num2str(t)]);

w_beamlets = ones(numBeamlet,1);
numBeam=numel(unique(beam_i_list));
for beam_i = 1:numBeam % assign weights to beamlets
    % beamlets from same beam get same initial weights
    w_beamlets(beam_i_list == beam_i) = w_beam(beam_i);
end
w_beamlets = w_beamlets + (2*rand(size(w_beamlets))-1) *0.1 .*w_beamlets; % small random perturbation
w_beamlets = 1.1* w_beamlets; % this just kicks the beamlets up a bit to make it easier for the optimizer to start

% -- GET FULL BEAMLET WEIGHTS --
options.MaxFunctionEvaluations = N_fcallback2;
% tic
fprintf('\n running full optimizer:')
w_fin = fmincon(fun2,w_beamlets,A,b,Aeq,beq,lb,ub,nonlcon,options);
fprintf('  done!:')
t=toc;
disp(['Optimization time for beamlets = ',num2str(t)]);


%% evaluate the results
D_full = reshape(beamlets * w_fin, size(Geometry.data));

%% save outputs
NLP_result.dose = D_full;
NLP_result.weights = w_fin;
save([patient_dir '\matlab_files\NLP_result.mat'], 'NLP_result');

plot_DVH(D_full, optGoal, optGoal_idx, targetMinMax_idx)
colorwash(Geometry.data, D_full, [-500, 500], [0, 110]);
% plot_DVH_robust(D_full, optGoal, optGoal_idx)
end

%% support functions
% ---- PENALTY FUNCTION ----
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 = sss
            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
% ------ supp: penalty for single scenario ------
function penalty = eval_f(x, optGoal, nrs_i, sss_i, rrs_i)
    penalty = 0;
    % for each condition
    for goal_i = 1:numel(optGoal)
        switch optGoal{goal_i}.function
            % min, max, min_sq, max_sq, LeastSquare, min_perc_Volume, max_perc_Volume
            case 'min'
                % penalize if achieved dose is lower than target dose
                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 = 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 'min_sq'
                % penalize if achieved dose is higher than target dose
                temp1=min(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));
                d_penalty = 1.0e0 * sum(temp1.^2);
            case 'max_sq'
                % penalize if achieved dose is higher than target dose
                temp1=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));
                d_penalty = 1.0e0 * sum(temp1.^2);
            case 'LeastSquare'
                % penalize with sum of squares any deviation from target
                % dose
                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);
            case 'min_perc_Volume'
                % penalize by amount of volume under threshold
                perc_vox = numel(find((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) > 0)) ...
                    / numel(optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.target);
                d_penalty = 3.0e5 * min(perc_vox-0.05, 0)
                
            case 'max_perc_Volume'
                % penalize by amount of volume under threshold
                perc_vox = numel(find((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) < 0)) ...
                    / numel(optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.target);
                d_penalty = 3.0e4 * min(perc_vox-0.05, 0)
                    
        end
        penalty = penalty + d_penalty * optGoal{goal_i}.opt_weight;
    end
end


% ---- MAKE ROI ROBUST ----
function optGoal = make_robust_optGoal(optGoal, RO_params, beamlets);
    % take regular optimal goal and translate it into several robust cases
    
    % RO_params - should have the information below
    % nrs - random scenarios
    % sss - system setup scenarios
    % rrs - random range scenarios
    
    % X - X>0 moves image right
    % Y - Y>0 moves image down
    % Z - in/out.
    
    shift_mag = 2; % vox of shift
    nrs_scene_list={[0,0,0]};

    
    % ----====#### CHANGE ROBUSTNESS HERE ####====----
    
%     sss_scene_list={[0,0,0]};
    sss_scene_list={[0,0,0], [-shift_mag,0,0], [shift_mag,0,0], [0,-shift_mag,0], [0,shift_mag,0], [0,0,-1], [0,0,1]};
%     sss_scene_list={[0,0,0], [-shift_mag,0,0], [shift_mag,0,0], [0,-shift_mag,0], [0,shift_mag,0],...
%         [-shift_mag*2,0,0], [shift_mag*2,0,0], [0,-shift_mag*2,0], [0,shift_mag*2,0]};

    % ----====#### CHANGE ROBUSTNESS HERE ####====----

    rrs_scene_list={[0,0,0]};

    [targetIn, meta] = nrrdread('C:\010-work\003_localGit\WiscPlan_v2\data\archive\CDP_data\CDP5_DP_target.nrrd');
%     [targetIn, meta] = nrrdread('C:\010-work\003_localGit\WiscPlan_v2\data\PD_HD_dicomPhantom\Tomo_DP_target.nrrd');
%     [targetIn, meta] = nrrdread('C:\010-work\003_localGit\WiscPlan_v2\data\archive\CDP_data\CDP5_DP_target.nrrd');
    
    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)
        % get target
        idx=optGoal{goal_i}.ROI_idx;
        targetImg1=zeros(optGoal{goal_i}.imgDim);
        targetImg1(idx)=1;
        % get beamlets
        
        for nrs_i = 1:optGoal{goal_i}.NbrRandScenarios          % num. of random scenarios
            % modify target and beamlets
            targetImg2=targetImg1;
            % beamlets stay the same
            
            for sss_i = 1 :optGoal{goal_i}.NbrSystSetUpScenarios   % syst. setup scenarios = sss
                % modify target and beamlets
                [targetImg3 idxValid]=get_RO_sss(targetImg2, sss_scene_list{sss_i});
                % beamlets stay the same
                
                for rrs_i = 1:optGoal{goal_i}.NbrRangeScenarios   % range scenario = rrs
                    % modify target and beamlets
                    targetImg4=targetImg3;
                    % beamlets stay the same
                    
                    %% make new target and beamlets
                    ROI_idx=[];
                    ROI_idx=find(targetImg4>0);
                    
                    if isfield(optGoal{goal_i}, 'target_alpha')
                        target = optGoal{goal_i}.target;
                        target = target(idxValid);
%                         disp('exists')
                    else
                        target = ones(numel(ROI_idx), 1) .* optGoal{goal_i}.D_final;
%                         disp('not exists')
                    end
                    
                    
                    beamlets_pruned = beamlets(ROI_idx, :);
                    
                    % save to optGoal output
                    optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.ROI_idx            = ROI_idx;
                    optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.beamlets_pruned    = beamlets_pruned;
                    optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.target             = target;
                end
            end
        end
    end
end
% ------ supp: RO case SSS ------
function [targetImg3 ia]=get_RO_sss(targetImg2, sss_scene_shift);
    % translate the target image 
    
    targetImg3 = imtranslate(targetImg2,sss_scene_shift);
    
    % now we need to figure out if any target voxels fell out during the
    % shift
    imgValid = imtranslate(targetImg3,-sss_scene_shift);
    imgInvalid = (targetImg2-imgValid);
    
    idx_1 = find(targetImg2);
    idx_2 = find(imgInvalid);
    
    [idxValid,ia] = setdiff(idx_1,idx_2);
    
    [C,ia, ib] = intersect(idx_1,idxValid);
    
end