function RDXTPS_optimSetup(Nbeamlets, patient_dir, Geometry)

% Sets up a dose prescription and creates the optimizer input files for
% proton plans.

if isempty(patient_dir)
    patient_dir = uifile('getdir', 'Select the patient directory');
end

if isempty(Geometry)
    load(fullfile(patient_dir, 'matlab_files', 'Geometry.mat'));
end

%% Setup dialog box
prompt = {...
    sprintf('Maximum number of iterations:'), ...
    sprintf('Save result every ? interations:'), ...
    sprintf('Modulation factor:\n(typically 5 for SS, 2 for DET)') ...
    };
defAns = {'2000', '200', '5'};

%% Prompt input dialog
answer = inputdlg(prompt, 'Optimization parameters', 1, defAns);

Niterations = str2double(answer{1});
Nperbatch = str2double(answer{2});
modFactor = str2double(answer{3});

%%
optSettings = [];

% set up the optimization: 
% All folder specified relative to optimizer executable
optSettings.optInfo.saveOptInfo = 'yes';
optSettings.optInfo.Niterations = Niterations;  % total number of iterations
optSettings.optInfo.Nperbatch = Nperbatch;     % number of iterations per batch
optSettings.optInfo.optFolder = patient_dir;
optSettings.optInfo.inputFile = 'optInput.txt'; % goes in the optimizer folder
optSettings.optInfo.inputFolder = 'opt_input';
optSettings.optInfo.outputFolder = 'opt_output';
optSettings.optInfo.modFactor = modFactor;
% Flag telling the optimizer whether to calculate an initial guess itself
% or to read in the initial guess file
optSettings.initialGuessInfo.calcInitialGuessFlag = 1; 

% DCW - the initial guess only works with beamlets stored in the matlab
% format as is done with the old matlab dose calculator.  thus for this
% application, initial guess will always be set to 'no'.
% initial guess calculation sometimes time consuming, turn off with flag: 
calculateInitialGuess = 'no';  

optSettings.prescInfo.savePresc = 'yes';

optSettings.beamletInfo.saveBeamletHeader = 'yes';
optSettings.beamletInfo.beamletFolder = 'beamlet_batch_files';
optSettings.beamletInfo.Nbeamlets = Nbeamlets;
optSettings.beamletInfo.Nbeamletbatches = 1;

%--------------------------------------------------------------------------
% DCW 10-13-07
% XMO 11-06-09

ROI_names = cellfun(@(c)c.name, Geometry.ROIS, 'UniformOutput', false);

[target_idx okay] = listdlg('ListString', ROI_names, ...
    'SelectionMode', 'single', 'Name', 'Target Selection', ...
    'PromptString', 'Please select the target ROI. ');
% if okay ~= 1
%     msgbox('Optimization setup aborted');
%     return;
% end

% parameters set to 1 if ROI is target, 0 otherwise
prescTable = cell(length(Geometry.ROIS),13);
for roi_ind = 1:length(Geometry.ROIS)
    prescTable{roi_ind,1} = Geometry.ROIS{roi_ind}.name;
    prescTable{roi_ind,2} = Geometry.ROIS{roi_ind}.name;
    % set all parameters to zero
    for iloop = 3:13
        prescTable{roi_ind,iloop} = 0;
    end
    % set default values for target
    if roi_ind == target_idx
        prescTable{roi_ind, 3} = 100;   % alpha
        prescTable{roi_ind, 4} = 1;     % max dose penalty
        prescTable{roi_ind, 5} = 60;    % max dose
        prescTable{roi_ind, 6} = 1;     % min dose penalty
        prescTable{roi_ind, 7} = 60;    % min dose
        prescTable{roi_ind,12} = 60;    % DVH under dose (prescription dose)
        prescTable{roi_ind,13} = 98;    % DVH under dose V (prescription volume)
    end
end

%--------------------------------------------------------------------------
% % load downsampled geometry file
% load(geometryFile);

% size of CT grid
[M,N,Q] = size(Geometry.rhomw);

% Create the prescription structure
presc = [];   % erase the old prescription structure, if it exists
presc.siz = [M N Q];  % dimensions of the prescription and tissue masks

ROIFlags = zeros(1,size(prescTable,1));  % flags to ensure all ROIs are found

% assign the parameters to the prescription structure
for j=1:size(prescTable,1)
    % pull out the prescription indices for each tissue
    tissName = prescTable{j,2};  % current tissue name

    % find the current tissue in the Geometry.ROIS cells
    for k=1:length(Geometry.ROIS)
        if strcmp(tissName,strrep(Geometry.ROIS{k}.name,'%',''))
            ind = Geometry.ROIS{k}.ind;  % extract ROI indices
            break;
        end
    end
    
    presc.tissue(j).name = prescTable{j,1};
    presc.tissue(j).ind = ind;
    presc.tissue(j).alpha = prescTable{j,3};
    presc.tissue(j).betaPlus = prescTable{j,4};
    presc.tissue(j).dPlus = single(zeros(size(ind)) + prescTable{j,5});
    presc.tissue(j).betaMinus = prescTable{j,6};
    presc.tissue(j).dMinus = single(zeros(size(ind)) + prescTable{j,7});
    presc.tissue(j).betaVPlus = prescTable{j,8};
    presc.tissue(j).dVPlus = prescTable{j,9};
    presc.tissue(j).vPlus = prescTable{j,10};
    presc.tissue(j).betaVMinus = prescTable{j,11};
    presc.tissue(j).dVMinus = prescTable{j,12};
    presc.tissue(j).vMinus = prescTable{j,13};
    presc.tissue(j).dPlus = single(zeros(size(ind)) + prescTable{j,5});
    presc.tissue(j).dMinus = single(zeros(size(ind)) + prescTable{j,7});
    ROIFlags(j) = 1;
end

% ensure that all of the ROIS were found
if sum(ROIFlags) ~= length(ROIFlags)
    error('Missed one or more ROIS.\n');
else
    fprintf('All ROIs located successfully for downsampled prescription\n');
end

optSettings.prescInfo.presc = presc;
optSettings.beamletInfo.beamletDim = [M N Q];

RDX_linlsqOptimization(optSettings);