function [optSettings, optResults] = loadOptResults(varargin)
% Loads the weights and dose distribution corresponding to a given input
% filename.  The input filename has the following format:
%bemal
% [S,R] = loadOptResults(inputFileName) loads all of the files from
%        the linlsq optimization associated with inputFileName.  DVHs are
%        calculated for all of the associated dose files.
% [S,R] = loadOptResults(inputFileName,'last') loads only the last
%        dose and beamlet weights files, and calculates the DVHs for each
%        tissue type only for those files.
%
% inputFileName has the following format:
% Niterations
% 2000
% Nperbatch
% 50
% prescription_filename
% input0/prescription.txt
% initial_beam_weights_filename
% input0/init_beam_weights.img
% beamlet_header_file
% beamletbatches0/beamlet_header.txt
% dose_batch_base_name
% output0/dosebatch
% dose_batch_extension
% img
% weight_batch_base_name
% output0/weightbatch
% weight_batch_extension
% img
% obj_func_name
% output0/objFunc.img
%
% Where all paths are relative to the directory containing the
% inputFileName argument.
%
% RTF 1/6/07

% constants
Ndvhbins = 1000; % number of bins to use for the cumulative DVH calculation

if length(varargin) == 1
    optInputFile = varargin{1};
    readAllOutputFiles = 1;
elseif length(varargin) == 2
    optInputFile = varargin{1};
    if strcmp(deblank(lower(varargin{2})),'last')
        readAllOutputFiles = 0;
    else
        readAllOutputFiles = 1;
    end
else
    error('Too many input arguments');
end

% Extract the folder name that contains the input file
inputFileNameRev = fliplr(optInputFile);  % flip the input filename around
% pop off the reversed file name
[fileNameRev,inputFolderRev] = strtok(inputFileNameRev,{'/','\'});
inputFolder = fliplr(inputFolderRev);

[optSettings,missingInfo] = loadOptSettings(optInputFile);  % load the optimization settings

if numel(missingInfo)
    fprintf('Information missing from optimization files:\n');
    for k=1:numel(missingInfo)
        fprintf('%s\n',missingInfo{k});
    end
    fprintf('\n');
end

% find which dose and weights files are present
Nbatches = ceil(optSettings.optInfo.Niterations/optSettings.optInfo.Nperbatch)+1;  % extra '1' is for initial guess
doseFileNames = cell(1,Nbatches);
doseFileExist = zeros(1,Nbatches);  % flags to test for existence of files
weightFileNames = cell(1,Nbatches);
weightFileExist = zeros(1,Nbatches);

for k=1:Nbatches
    doseFileNames{k} = [inputFolder optSettings.optInfo.outputFolder '/' ...
        optSettings.optInfo.doseBatchBaseName ...
        num2str((k-1)*optSettings.optInfo.Nperbatch) '.' ...
        optSettings.optInfo.doseBatchExtension];
    weightFileNames{k} = [inputFolder optSettings.optInfo.outputFolder '/' ...
        optSettings.optInfo.weightBatchBaseName ...
        num2str((k-1)*optSettings.optInfo.Nperbatch) '.' ...
        optSettings.optInfo.weightBatchExtension];

    fid = fopen(doseFileNames{k},'rb');
    if fid == -1  % batch file doesn't exist
        doseFileExist{k} = 0;
    else
        fclose(fid);
        doseFileExist(k) = 1;  % mark that file exists
    end

    fid = fopen(weightFileNames{k},'rb');
    if fid == -1  % batch file doesn't exist
        weightFileExist{k} = 0;
    else
        fclose(fid);
        weightFileExist(k) = 1;  % mark that file exists
    end
end

optResults = [];

if readAllOutputFiles == 1
    optResults.dose = {};
    optResults.weights = {};
    % read in the appropriate dose files
    for k=1:length(doseFileNames)
        if doseFileExist(k) == 1
            fid = fopen(doseFileNames{k},'rb');
            dose = fread(fid,'float');
            fclose(fid);
            dose = reshape(dose,optSettings.prescInfo.presc.siz);
            optResults.dose{k} = dose;
        end
    end

    % read in the appropriate weights files
    for k=1:length(weightFileNames)
        if weightFileExist(k) == 1
            fid = fopen(weightFileNames{k},'rb');
            weights = fread(fid,'float');
            fclose(fid);
            optResults.weights{k} = weights;
        end
    end
else  % read in only the last dose and weights files
    for k=length(doseFileNames):-1:0
        if doseFileExist(k) == 1  % found last dose file
            break;
        end
    end
    
    fid = fopen(doseFileNames{k},'rb');
    dose = fread(fid,'float');
    fclose(fid);
    dose = reshape(dose,optSettings.prescInfo.presc.siz);
    optResults.dose = dose;
    
    fid = fopen(weightFileNames{k},'rb');
    weights = fread(fid,'float');
    fclose(fid);
    optResults.weights = weights;
end

optResults.presc = optSettings.prescInfo.presc;

% calculate a max dose vector
if readAllOutputFiles == 1
    for k=1:length(optResults.dose)
        dmax(k) = 1.1*max(optResults.dose{k}(:));
    end
    % calculate DVHs for each tissue, for each batch number
    for m=1:length(optResults.presc.tissue)
        tissMask = zeros(optResults.presc.siz,'int8');
        tissMask(optResults.presc.tissue(m).ind) = 1;
        for k=1:length(optResults.dose)
            dvhbins = [0:Ndvhbins-1]*dmax(k)/Ndvhbins;
            optResults.presc.tissue(m).dvhbins{k} = dvhbins;
            optResults.presc.tissue(m).dvh{k} ...
                = dvh(optResults.dose{k},tissMask,dvhbins)';
        end
    end  
else
    dmax = 1.1*max(optResults.dose(:));
    % calculate DVHs for each tissue
    for m=1:length(optResults.presc.tissue)
        tissMask = zeros(optResults.presc.siz,'int8');
        tissMask(optResults.presc.tissue(m).ind) = 1;
        dvhbins = [0:Ndvhbins-1]*dmax/Ndvhbins;
        optResults.presc.tissue(m).dvhbins = dvhbins;
        optResults.presc.tissue(m).dvh ...
            = dvh(optResults.dose,tissMask,dvhbins)';
    end
end

% load the objective function
fid = fopen([inputFolder '/' optSettings.optInfo.outputFolder '/' optSettings.optInfo.objFuncFileName],'rb');
optResults.objFunc = fread(fid,'float');
fclose(fid);