WiscPlan/WiscPlanPhotonkV125/matlab_frontend/lab/RDX_linlsqOptimization.m

function [weights, dose] = RDX_linlsqOptimization(optSettings)
% Input:  optSettings structure, defined as follows:
% optSettings
%  .prescInfo
%      .savePresc = 'yes' or 'no' (defaults to 'no')
%      .presc (prescription structure, which has the form:
%              .presc.siz = [M N Q] size of each beamlet
%              .presc.tissue = Structure array of tissue-specific
%                  prescription.  savePrescription subfunction has details.
%              Defaults to an empty matrix. 
%      .prescFileName (defaults to 'prescription.txt')
% 
%  .beamletInfo
%      .saveBeamlets = 'yes' or 'no' (defaults to 'no')
%      .saveBeamletHeader = 'yes' or 'no' (defaults to 'no')
%      .beamletFolder (defaults to 'beamlets')
%      .beamlets (cell array of Matlab-style sparse beamlets stored as:
%                 .beamlets{n}.Dij{m}, where Dij contains the sparse data)
%                 Defaults to an empty matrix.
%      .beamletBatchBaseName  (defaults to 'beamletbatch')
%      .beamletBatchExtension (defaults to 'bin')
%      .beamletHeaderFileName (defaults to 'beamlet_header.txt')
%      .Nbeamlets             (defaults to 1)
%      .Nbeamletbatches       (defaults to 1)
%      .beamletDim            (defaults to [1 1 1])
%
%  .optInfo
%      .saveOptInfo = 'yes' or 'no' (defaults to 'no')
%      .optRun = 'yes' or 'no'  (defaults to 'no')
%      .waitForResults = 'yes' or 'no'  (defaults to 'no')
%      .optFolder (Folder containing optimizer's Visual C++ code)
%                 Defaults to 'C:\rtflynn\library\matlabAdditions\optimization\linlsq\linlsqOpt'
%      .inputFolder (defaults to 'input')
%      .outputFolder (defaults to 'output')
%      .inputFile (defaults to optInput.txt)
%      .Niterations (Defaults to 100)
%      .Nperbatch   (Defaults to 50)
%      .doseBatchBaseName (defaults to 'dosebatch')
%      .doseBatchExtension (defaults to 'img')
%      .weightBatchBaseName (defaults to 'weightbatch')
%      .weightBatchExtension (defaults to 'img')
%      .optExeName (defaults to 'linlsqOpt.exe')
%      .objFuncFileName (defaults to 'objFunc.bin')
%      .modFactor (defaults to 1e10)
% 
%  .initialGuessInfo
%      .saveInitialGuess = 'yes' or 'no' (defaults to 'no')
%      .w0 (initial guess for all beamlet weights, defaults to empty matrix)
%      .initBeamWeightFile (defaults to 'init_beam_weight_file')
%      .calcInitialGuessFlag (defaults to 1)  % if 1, optimizer calculates
%            the initial guess
%
% Outputs:
% w              optimized beamlet weights
% dose           optimized dose distribution
%
% Outputs are empty matrices unless the 'waitForResults' flag is set to
% 'yes'.  Otherwise, load results at a later time.
%
% Ryan Flynn, David Westerly
% Xiaohu Mo
% Fixed pathsep (/\) portability issue

% Check input structure for inconistencies, set defaults for fields that
% were optional and not included.
optSettings = parseInput(optSettings);

% From this point on, it is assumed that all fields that are required for
% any functions for the rest of this script are present.

origDir = pwd;  % copy the directory path from which this script was called

warning off;

% Change to the directory containing the linlsqOpt source code.
cd(optSettings.optInfo.optFolder);

% save prescription info
savePrescription(optSettings);

% save beamlet info
saveBeamlets(optSettings);

% save initial guess info
saveInitialGuess(optSettings);

% save optimization info
saveOptInfo(optSettings);

% run optimizer
[weights, dose] = runOptimizer(optSettings);

cd(origDir);   % go back to the original directory


function optSettings = parseInput(optSettings)
% Checks the optSettings input structure for inconsistencies and applies
% defaults to optional fields that were not included in the input
% structure.  

% set optimizer defaults
optSettingsDefault.prescInfo.savePresc = 'no';
optSettingsDefault.prescInfo.presc = [];
optSettingsDefault.prescInfo.prescFileName = 'prescription.txt';
optSettingsDefault.beamletInfo.saveBeamlets = 'no';
optSettingsDefault.beamletInfo.saveBeamletHeader = 'no';
optSettingsDefault.beamletInfo.beamletFolder = 'beamlets';
optSettingsDefault.beamletInfo.beamlets = [];
optSettingsDefault.beamletInfo.beamletBatchBaseName = 'beamletbatch';
optSettingsDefault.beamletInfo.beamletBatchExtension = 'bin';
optSettingsDefault.beamletInfo.beamletHeaderFileName = 'beamlet_header.txt';
optSettingsDefault.beamletInfo.Nbeamlets = 1;
optSettingsDefault.beamletInfo.Nbeamletbatches = 1;
optSettingsDefault.beamletInfo.beamletDim = [1 1 1];
optSettingsDefault.optInfo.saveOptInfo = 'no';
optSettingsDefault.optInfo.optRun = 'no';
optSettingsDefault.optInfo.waitForResults = 'no';
optSettingsDefault.optInfo.optFolder = ...
    'F:\Treatment_Planning\optimizer';
optSettingsDefault.optInfo.inputFolder = 'input';
optSettingsDefault.optInfo.outputFolder = 'output';
optSettingsDefault.optInfo.inputFile = 'optInput.txt';
optSettingsDefault.optInfo.Niterations = 100;
optSettingsDefault.optInfo.Nperbatch = 50;
optSettingsDefault.optInfo.doseBatchBaseName = 'dosebatch';
optSettingsDefault.optInfo.doseBatchExtension = 'img';
optSettingsDefault.optInfo.weightBatchBaseName = 'weightbatch';
optSettingsDefault.optInfo.weightBatchExtension = 'img';
optSettingsDefault.optInfo.optExeName = 'linlsqOptimizer.exe';
optSettingsDefault.optInfo.objFuncFileName = 'objFunc.bin';
optSettingsDefault.optInfo.modFactor = 1e10;  % modulation factor
optSettingsDefault.initialGuessInfo.saveInitialGuess = 'no';
optSettingsDefault.initialGuessInfo.w0 = [];
optSettingsDefault.initialGuessInfo.initBeamWeightFile = 'init_beam_weight_file.img';
optSettingsDefault.initialGuessInfo.calcInitialGuessFlag = 1;

% Traverse the optSettings structure tree, and, whenever a leaf is missing
% in the tree, fill it in with its default value.
nodeVec = [1 1];   % vector corresponding to nodes in the structure tree
while length(nodeVec) > 1  % traverse as long as we're not back at the top of the tree
    % Create a field access string, which is used to access the node
    % referred to by nodeVec.
    fldAccStr = 'optSettingsDefault';   
    for k=2:length(nodeVec)  
        eval(['fieldNames = fieldnames(' fldAccStr ');']);
        fldAccStr = [fldAccStr '.' fieldNames{nodeVec(k)}];
    end
    
    % Check if the field is a number or a string (leaf).
    eval(['numFlag = isnumeric(' fldAccStr ');']);
    eval(['charFlag = ischar(' fldAccStr ');']);
    
    if numFlag == 1 | charFlag == 1  % found a leaf
        % Traverse input structure, searching for leaves that exist in the
        % default structure but do not exist in the input structure.  This
        % is done by working down the input tree until either a node or a
        % leaf is missing from the input
        fldAccStrInp = 'optSettings';
        fldAccStrDummy = 'optSettingsDefault';
        for k=2:length(nodeVec)
            eval(['fieldNamesDummy = fieldnames(' fldAccStrDummy ');']);
            eval(['fieldFlag = isfield(' fldAccStrInp ',''' fieldNamesDummy{nodeVec(k)} ''');']);
            if fieldFlag == 1 
                fldAccStrInp = [fldAccStrInp '.' fieldNamesDummy{nodeVec(k)}];
                fldAccStrDummy = [fldAccStrDummy '.' fieldNamesDummy{nodeVec(k)}];
            else  % Leaf cannot exist, so set it to its default.
                [tok,rem] = strtok(fldAccStr,'.');
                fldAccStrInp = ['optSettings' rem];
                eval([fldAccStrInp ' = ' fldAccStr ';']);
                break;
            end
        end
        
        % go to the next node vector
        while length(nodeVec) > 1
            fldAccStr = 'optSettingsDefault';
            for k=2:length(nodeVec)
                eval(['fieldNames = fieldnames(' fldAccStr ');']);
                fldAccStr = [fldAccStr '.' fieldNames{nodeVec(k)}];
            end
            % Go to next leaf or nodes on the same level, if one exists.
            if nodeVec(end) + 1 <= length(fieldNames)
                nodeVec(end) = nodeVec(end) + 1;
                break;
            else % No more leaves or nodes on this level, so go up a level.
                nodeVec(:,end) = []; 
            end
        end
    else  % found another field that is not a leaf
        nodeVec = [nodeVec 1];  % go down a level in the tree
    end
end

% If beamlets were included with optSettings, ensure that the
% Nbeamletbatches, Nbeamlets, and beamletDim fields are consistent
if ~isempty(optSettings.beamletInfo.beamlets)
    % get the beamlet dimensions
    [Xcount,Ycount,Zcount] = size(full(optSettings.beamletInfo.beamlets{1}.Dij{1}));
    optSettings.beamletInfo.beamletDim = [Xcount Ycount Zcount];
    optSettings.beamletInfo.Nbeamletbatches = numel(optSettings.beamletInfo.beamlets);
    optSettings.beamletInfo.Nbeamlets = 0;
    for k=1:optSettings.beamletInfo.Nbeamletbatches
        optSettings.beamletInfo.Nbeamlets = optSettings.beamletInfo.Nbeamlets ...
            + numel(optSettings.beamletInfo.beamlets{k}.Dij);
    end
else
    % If optSettings.prescInfo.presc.siz was defined, then that field
    % should contain beamlet dimensions.
    if isfield(optSettings.prescInfo.presc,'siz')
        optSettings.beamletInfo.beamletDim = optSettings.prescInfo.presc.siz;
    end
end

function savePrescription(optSettings)
% If the savePresc field is 'yes', saves a prescription structure to a file 
% format that is useable by the linear least squares optimizer.
%% .prescInfo
%%     .savePresc = 'yes' or 'no'
%      .presc (prescription structure, which has the form:
%              .presc.siz = [M N Q] size of each beamlet
%              .presc.tissue = structure array of tissue-specific
%                  prescription.  savePrescription subfunction has details.
%      .prescFile (defaults to 'prescription.txt')

infValue = 10e10;   % value of infinity

% Check the savePresc flag.  If 'no', do not save prescription
if strcmp(lower(optSettings.prescInfo.savePresc),'no')
    return;
end

% create input and output folders
success = mkdir(optSettings.optInfo.inputFolder); 

% delete old output folder if it exists first
success = rmdir(optSettings.optInfo.outputFolder,'s');
success = mkdir(optSettings.optInfo.outputFolder); 

% Check the prescription structure to ensure all proper fields are present.
prescTissueFields = {'name' 'alpha' 'betaVPlus' 'dVPlus' 'vPlus' 'betaVMinus' ...
            'dVMinus' 'vMinus' 'betaPlus' 'betaMinus' 'ind' 'dPlus' 'dMinus'};
            
if isempty(optSettings.prescInfo.presc)
    error('optSettings.prescInfo.presc not defined -- define it or set optSettings.prescInfo.savePresc to ''no''');
end
            
if ~isfield(optSettings.prescInfo.presc,'siz')
    error('Missing ''siz'' field from optSettings.prescInfo.presc structure');
end

% siz is present, so make sure it is a 3 element vector
if numel(optSettings.prescInfo.presc.siz) ~= 3
    error('The ''siz'' field from optSettings.prescInfo.presc structure must be a 3-element vector');
end

if ~isfield(optSettings.prescInfo.presc,'tissue')
    error('Missing ''tissue'' field from prescInfo.presc structure');
end

% search for subfields of tissue field
for k=1:numel(prescTissueFields)
    if ~isfield(optSettings.prescInfo.presc.tissue,prescTissueFields{k})
        error(['Missing ''' prescTissueFields{k} ''' field from optSettings.prescInfo.presc.tissue structure']); 
    end
end

% write the prescription file
fid = fopen([optSettings.optInfo.inputFolder filesep optSettings.prescInfo.prescFileName],'w');

if fid == -1
    error('Failed to open prescription file.');
end

fprintf(fid,'Ntissue\n');
fprintf(fid,[num2str(numel(optSettings.prescInfo.presc.tissue)) '\n\n']);

for k=1:numel(optSettings.prescInfo.presc.tissue)   % write prescription parameters for each tissue type
    % write the current tissue number
    fprintf(fid,'tissueNum\n');
    fprintf(fid,[num2str(k-1) '\n']);
    fprintf(fid,'name\n');
    fprintf(fid,[optSettings.prescInfo.presc.tissue(k).name '\n']);
    
    % throughout writing process, always check for infinities before printing values
    
    % write the tissue importance
    fprintf(fid,'alpha\n');
    if isinf(optSettings.prescInfo.presc.tissue(k).alpha)
        fprintf(fid,[num2str(infValue) '\n']);
    else
        fprintf(fid,[num2str(optSettings.prescInfo.presc.tissue(k).alpha) '\n']);
    end

    fprintf(fid,'betaVPlus\n');
    if isinf(optSettings.prescInfo.presc.tissue(k).betaVPlus)
        fprintf(fid,[num2str(infValue) '\n']);
    else
        fprintf(fid,[num2str(optSettings.prescInfo.presc.tissue(k).betaVPlus) '\n']);
    end
    
    fprintf(fid,'dVPlus\n');
    if isinf(optSettings.prescInfo.presc.tissue(k).dVPlus)
        fprintf(fid,[num2str(infValue) '\n']);
    else
        fprintf(fid,[num2str(optSettings.prescInfo.presc.tissue(k).dVPlus) '\n']);
    end
    
    fprintf(fid,'vPlus\n');
    if isinf(optSettings.prescInfo.presc.tissue(k).vPlus)
        fprintf(fid,[num2str(infValue) '\n']);
    else
        fprintf(fid,[num2str(optSettings.prescInfo.presc.tissue(k).vPlus) '\n']);
    end
    
    fprintf(fid,'betaVMinus\n');
    if isinf(optSettings.prescInfo.presc.tissue(k).betaVMinus)
        fprintf(fid,[num2str(infValue) '\n']);
    else
        fprintf(fid,[num2str(optSettings.prescInfo.presc.tissue(k).betaVMinus) '\n']);
    end
    
    fprintf(fid,'dVMinus\n');
    if isinf(optSettings.prescInfo.presc.tissue(k).dVMinus)
        fprintf(fid,[num2str(infValue) '\n']);
    else
        fprintf(fid,[num2str(optSettings.prescInfo.presc.tissue(k).dVMinus) '\n']);
    end
    
    fprintf(fid,'vMinus\n');
    if isinf(optSettings.prescInfo.presc.tissue(k).vMinus)
        fprintf(fid,[num2str(infValue) '\n']);
    else
        fprintf(fid,[num2str(optSettings.prescInfo.presc.tissue(k).vMinus) '\n']);
    end
    
    
    fprintf(fid,'betaPlus\n');
    if isinf(optSettings.prescInfo.presc.tissue(k).betaPlus)
        fprintf(fid,[num2str(infValue) '\n']);
    else
        fprintf(fid,[num2str(optSettings.prescInfo.presc.tissue(k).betaPlus) '\n']);
    end
    
    % create a name for the underdose penalty prescription file
    dose_presc_filename = [regexprep(optSettings.prescInfo.presc.tissue(k).name,{' ','\','/'},'') 'DosePlus.bin'];
    fprintf(fid,'dosePlusFilename\n');
    fprintf(fid,'%s\n', fullfile(optSettings.optInfo.optFolder, optSettings.optInfo.inputFolder, dose_presc_filename));
    
    % save the mask as a sparse matrix
    fid2 = fopen([optSettings.optInfo.inputFolder filesep dose_presc_filename],'wb');
    fwrite(fid2,optSettings.prescInfo.presc.siz,'int');
    fwrite(fid2,length(optSettings.prescInfo.presc.tissue(k).ind),'int');
    fwrite(fid2,optSettings.prescInfo.presc.tissue(k).ind-1,'int');
    fwrite(fid2,optSettings.prescInfo.presc.tissue(k).dPlus,'single');
    fclose(fid2);
    
    fprintf(fid,'betaMinus\n');
    if isinf(optSettings.prescInfo.presc.tissue(k).betaMinus)
        fprintf(fid,[num2str(infValue) '\n']);
    else
        fprintf(fid,[num2str(optSettings.prescInfo.presc.tissue(k).betaMinus) '\n']);
    end
    
    % create a name for the overdose penalty prescription file
    dose_presc_filename = [regexprep(optSettings.prescInfo.presc.tissue(k).name,{' ','\','/'},'') 'DoseMinus.bin'];
    fprintf(fid,'doseMinusFilename\n');
    fprintf(fid,'%s\n', fullfile(optSettings.optInfo.optFolder, optSettings.optInfo.inputFolder, dose_presc_filename));
    
    % save the mask as a sparse matrix
    fid2 = fopen([optSettings.optInfo.inputFolder filesep dose_presc_filename],'wb');
    fwrite(fid2,optSettings.prescInfo.presc.siz,'int');
    fwrite(fid2,length(optSettings.prescInfo.presc.tissue(k).ind),'int');
    fwrite(fid2,optSettings.prescInfo.presc.tissue(k).ind-1,'int');
    fwrite(fid2,optSettings.prescInfo.presc.tissue(k).dMinus,'single');
    fclose(fid2);
    
    fprintf(fid,'\n');
end

fclose(fid);

function saveBeamlets(optSettings)
% Saves a set of beamlets to filename.  Dij is a cell array of beamlets,
% each of which must be sparse 2D Matlab matrices.
%
%% .beamletInfo
%%     .saveBeamlets = 'yes' or 'no' (if 'yes', beamlets saved)
%%     .beamletFolder (required, folder containing beamlets their header)
%      .beamlets (cell array of Matlab-style sparse beamlets stored as:
%                 .beamlets{n}.Dij{m}, where Dij contains the sparse data)
%      .beamletBatchBaseName (defaults to 'beamletbatch')
%      .beamletBatchExtension (defaults to 'bin')
%      .beamletHeaderFileName (defaults to 'beamlet_header.txt')
%      .Nbeamlets
%      .Nbeamletbatches
%      .beamletDim

% Check if the beamlet header will need to be saved
if strcmp(lower(optSettings.beamletInfo.saveBeamletHeader),'no') ...
    && strcmp(lower(optSettings.beamletInfo.saveBeamlets),'no')
    return;
end

if strcmp(lower(optSettings.beamletInfo.saveBeamlets),'yes')
    % delete all of the old beamlets and header info
    success = rmdir([optSettings.optInfo.optFolder '/' optSettings.beamletInfo.beamletFolder],'s');
    mkdir([optSettings.optInfo.optFolder '/' optSettings.beamletInfo.beamletFolder]);
end

% create a folder for the beamlets if it doesn't already exist
mkdir(optSettings.beamletInfo.beamletFolder);

% create the beamlet header file
fid = fopen([optSettings.beamletInfo.beamletFolder '/' optSettings.beamletInfo.beamletHeaderFileName],'w');
if fid == -1
    error('Failed to open beamlet header file.');
end

fprintf(fid,'beamletDim\n');
fprintf(fid,'%g %g %g\n',optSettings.beamletInfo.beamletDim);
fprintf(fid,'Nbeamlets\n');
fprintf(fid,'%g\n',optSettings.beamletInfo.Nbeamlets);
fprintf(fid,'Nbeamletbatches\n');
fprintf(fid,'%g\n',optSettings.beamletInfo.Nbeamletbatches);
fprintf(fid,'beamletFolder\n');
fprintf(fid,'%s\n', fullfile(optSettings.optInfo.optFolder, optSettings.beamletInfo.beamletFolder));
fprintf(fid,'beamletBatchBaseName\n');
fprintf(fid,[optSettings.beamletInfo.beamletBatchBaseName '\n']);
fprintf(fid,'beamletBatchExtension\n');
fprintf(fid,[optSettings.beamletInfo.beamletBatchExtension '\n']);
fclose(fid);

% Check if the beamlets need to be saved
if strcmp(lower(optSettings.beamletInfo.saveBeamlets),'no')
    return;
end

if isempty(optSettings.beamletInfo.beamlets)
    error('No beamlets to save -- supply beamlets or set saveBeamlets to ''no''');
end

% Save the beamlets
beamletNum = 0;
for k=1:optSettings.beamletInfo.Nbeamletbatches
    beamletFileName = [optSettings.beamletInfo.beamletFolder '/' ...
        optSettings.beamletInfo.beamletBatchBaseName num2str(k-1) '.' ...
        optSettings.beamletInfo.beamletBatchExtension];
    fid = fopen(beamletFileName,'wb');
    if fid == -1
        error(['Failed to open ' beamletFileName]);
    end
    Nbib = numel(optSettings.beamletInfo.beamlets{k}.Dij);  % number of beamlet in batch
    fwrite(fid,Nbib,'int');
    for m=1:Nbib
        % dimensions of current beamlet
        [Nx,Ny,Nz] = size(full(optSettings.beamletInfo.beamlets{k}.Dij{m}));
        if Nx - optSettings.beamletInfo.beamletDim(1) ~= 0 ...
        || Ny - optSettings.beamletInfo.beamletDim(2) ~= 0 ...
        || Nz - optSettings.beamletInfo.beamletDim(3) ~= 0
            error(['Dimensions of optSettings.beamletInfo.beamlets{' ...
                num2str(k) '}.Dij{' num2str(m) '} are not (' ...
                num2str(optSettings.beamletInfo.beamletDim(1)) ',' ...
                num2str(optSettings.beamletInfo.beamletDim(2)) ',' ...
                num2str(optSettings.beamletInfo.beamletDim(3)) ')']);
        end
        % save the beamlet
        fwrite(fid,beamletNum,'int');  % beamlet number
        fwrite(fid,optSettings.beamletInfo.beamletDim,'int');  % beamlet dimensions
        ind = find(optSettings.beamletInfo.beamlets{k}.Dij{m}); % beamlet indices
        ind = ind - 1;  % convert Matlab indices to C indices
        fwrite(fid,numel(ind),'int');  % number of non-zero indices
        fwrite(fid,ind,'int');         % non-zero indices
        fwrite(fid,nonzeros(optSettings.beamletInfo.beamlets{k}.Dij{m}),'float');
        beamletNum = beamletNum + 1;
    end
    fclose(fid);  % close the current beamlet batch
end

function saveInitialGuess(optSettings)
% save the initial guess

% exit if the initial guess is not to be saved
if strcmp(lower(optSettings.initialGuessInfo.saveInitialGuess),'no')
    return;
end

% create input and output folders
success = mkdir(optSettings.optInfo.inputFolder); 
success = mkdir(optSettings.optInfo.outputFolder); 

initBeamWeightFile = [optSettings.optInfo.optFolder filesep optSettings.optInfo.inputFolder ...
    filesep optSettings.initialGuessInfo.initBeamWeightFile];

fid = fopen(initBeamWeightFile,'wb');
if fid == -1
    error(['Failed to open initial beam weight file: ' initBeamWeightFile]);
end
fwrite(fid,optSettings.initialGuessInfo.w0,'float');
fclose(fid);

function saveOptInfo(optSettings)

% Exit if optimization info is not supposed to be saved
if strcmp(lower(optSettings.optInfo.saveOptInfo),'no')
    return;
end

% create input and output folders
success = mkdir(optSettings.optInfo.inputFolder); 
success = mkdir(optSettings.optInfo.outputFolder); 

% create the optimization input file
optInfoFileName = optSettings.optInfo.inputFile;
fid = fopen(optInfoFileName,'w');
if fid == -1
    error(['Failed to open optimization input file: ' optInfoFileName]);
end

fprintf(fid,'Niterations\n');
fprintf(fid,'%d\n',optSettings.optInfo.Niterations);
fprintf(fid,'Nperbatch\n');
fprintf(fid,'%d\n',optSettings.optInfo.Nperbatch);
fprintf(fid,'prescFile\n');
fprintf(fid,'%s\n',fullfile(optSettings.optInfo.optFolder, optSettings.optInfo.inputFolder, optSettings.prescInfo.prescFileName));
fprintf(fid,'initBeamWeightFile\n');
if optSettings.initialGuessInfo.calcInitialGuessFlag == 0
    % supply the initial guess file
    fprintf(fid,'%s\n',[optSettings.optInfo.optFolder '/' optSettings.initialGuessInfo.initBeamWeightFile]);
else % no need to supply file, as optimizer will calculate initial guess
    fprintf(fid,'%s\n','optimizerCalculatesInitialGuess');
end
fprintf(fid,'beamletHeaderFileName\n');
fprintf(fid,'%s\n', fullfile(optSettings.optInfo.optFolder, optSettings.beamletInfo.beamletFolder, optSettings.beamletInfo.beamletHeaderFileName));
fprintf(fid,'doseBatchBaseName\n');
fprintf(fid,'%s\n', fullfile(optSettings.optInfo.optFolder, optSettings.optInfo.outputFolder, optSettings.optInfo.doseBatchBaseName));
fprintf(fid,'doseBatchExtension\n');
fprintf(fid,[optSettings.optInfo.doseBatchExtension '\n']);
fprintf(fid,'weightBatchBaseName\n');
fprintf(fid,'%s\n', fullfile(optSettings.optInfo.optFolder, optSettings.optInfo.outputFolder, optSettings.optInfo.weightBatchBaseName));
fprintf(fid,'weightBatchExtension\n');
fprintf(fid,[optSettings.optInfo.weightBatchExtension '\n']);
fprintf(fid,'objFuncFileName\n');
fprintf(fid,'%s\n', fullfile(optSettings.optInfo.optFolder, optSettings.optInfo.outputFolder, optSettings.optInfo.objFuncFileName));
fprintf(fid,'modFactor\n');
fprintf(fid,[num2str(optSettings.optInfo.modFactor) '\n']);
fprintf(fid,'calcInitialGuessFlag\n');
fprintf(fid,[num2str(optSettings.initialGuessInfo.calcInitialGuessFlag) '\n']);
fclose(fid);

function [weights, dose] = runOptimizer(optSettings)

dose = [];
weights = [];

% Exit function of optimizer is not to be run.
if strcmp(lower(optSettings.optInfo.optRun),'no')
    return;
end

currDir = pwd;  % save a copy of the current directory path

% copy the most recent build of optimizer to the optFolder
exeDebugFile = [optSettings.optInfo.optFolder '\Debug\' optSettings.optInfo.optExeName];
exeFile = [optSettings.optInfo.optFolder filesep optSettings.optInfo.optExeName];
success = copyfile(exeDebugFile,exeFile);

if success == 0
    error(['Error copying ' exeDebugFile ' to ' exeFile]);
end

% do the optimization
delete('optDone.txt');   % delete the old dose file if it exists
fprintf('Running optimizer...\n');
% run the optimization from the console
eval(['! ' optSettings.optInfo.optExeName ' ' optSettings.optInfo.inputFile ' &']); 

% Exit if the user did not wish to wait for results
if strcmp(lower(optSettings.optInfo.waitForResults),'no')
    return;
end

% Check for the flag indicating that the optimization is finished
while ~exist('optDone.txt','file')
    pause(1);   % pause before resuming
end

% Optimization finished, so load the beamlet weights
fid = fopen([optSettings.optInfo.outputFolder '/' ...
    optSettings.optInfo.weightBatchBaseName num2str(optSettings.optInfo.Niterations) '.' ...
    optSettings.optInfo.weightBatchExtension],'rb');
if fid == -1
    error('Unable to open optimized weight results.');
end

weights = fread(fid,'float');
fclose(fid);

% load the dose distribution
fid = fopen([optSettings.optInfo.outputFolder '/' ...
    optSettings.optInfo.doseBatchBaseName num2str(optSettings.optInfo.Niterations) '.' ...
    optSettings.optInfo.doseBatchExtension],'rb');
if fid == -1
    error('Unable to open optimized dose results.');
end

dose = fread(fid,'float');
try
    dose = reshape(dose,optSettings.beamletInfo.beamletDim);
catch
    error(['Unable to reshape optimized dose distribution. Ensure ' ...
     'that optSettings.beamletInfo.beamletDim field matches dose grid dimensions']);
end