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(optimizerFolder,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(optimizerFolder,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

warning off;

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

if length(varargin) == 1
    error('Too few input arguments');
elseif length(varargin) == 2
    optimizerFolder = varargin{1};
    optInputFile = varargin{2};
    readAllOutputFiles = 1;
elseif length(varargin) == 3
    optimizerFolder = varargin{1};
    optInputFile = varargin{2};
    if strcmp(deblank(lower(varargin{3})),'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);
if isempty(inputFolder)  % input folder is the current folder
    inputFolder = [pwd '/'];
end

% get the patient header
headerRev = fliplr(optimizerFolder);
% pop off patient name
[headerNameRev,optFolderRev] = strtok(headerRev,{'/','\'});
patientHeader = fliplr(headerNameRev);

[optSettings,missingInfo] = loadOptSettings(optimizerFolder,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 the dose grid size
if isfield(optSettings,'beamletInfo') && isfield(optSettings.beamletInfo,'beamletDim')
    siz = optSettings.beamletInfo.beamletDim;
elseif isfield(optSettings,'prescInfo') && isfield(optSettings.prescInfo,'presc') ...
        && isfield(optSettings.prescInfo.presc,'siz')
    siz = optSettings.prescInfo.presc.siz;
else
    siz = [];
end

if isfield(optSettings,'optInfo') && isfield(optSettings.optInfo,'Niterations') ...
        && isfield(optSettings.optInfo,'Nperbatch') && isfield(optSettings.optInfo,'outputFolder')...
        && isfield(optSettings.optInfo,'doseBatchBaseName') && isfield(optSettings.optInfo,'doseBatchExtension')...
        && isfield(optSettings.optInfo,'weightBatchBaseName') && isfield(optSettings.optInfo,'weightBatchExtension')
    % find which dose and weights files are present
    % must 'ls' the opt_output folder.  One file contains the objective
    % function, while the others are either dose batches or weight batch
    % files.

    % xmo: MUST use mod verision of ls in *nix
    atmp = lsl(optSettings.optInfo.outputFolder);
    btmp = size(atmp);
    if btmp(1)>4
        Nbatches = ceil((btmp(1) - 3)/2);
    elseif btmp(1)==4
        Nbatches = 1;
    else
        error('No dose batches found');
    end
    % 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);

    % Dave's new implementation
    k=0;
    for t=1:btmp(1)
        if ~isempty(findstr(optSettings.optInfo.doseBatchBaseName,mat2str(atmp(t,:))));
            k=k+1;
            batchname = strtrim(strrep(mat2str(atmp(t,:)),'''',''));
            str = strrep(strtok(batchname,'.'),optSettings.optInfo.doseBatchBaseName,'');
            batchNames{k} = str;
            doseFileNames{k} = [optSettings.optInfo.outputFolder '/' ...
                strtrim(strrep(mat2str(atmp(t,:)),'''','')) ];
            weightFileNames{k} = strrep(doseFileNames{k},optSettings.optInfo.doseBatchBaseName,optSettings.optInfo.weightBatchBaseName);
        end
    end

    for k=1:Nbatches
        % ryans implementation
        %         doseFileNames{k} = [optSettings.optInfo.outputFolder '/' ...
        %             optSettings.optInfo.doseBatchBaseName ...
        %             num2str((k-1)*optSettings.optInfo.Nperbatch) '.' ...
        %             optSettings.optInfo.doseBatchExtension];
        %         weightFileNames{k} = [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 = [];

    % assign names to optResults structure array
    optResults.batchNames = batchNames;

    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=>float');
                fclose(fid);
                if numel(siz) == 3 % reshape dose if siz definedre
                    dose = reshape(dose,siz);
                end
                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=>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=>float');
        fclose(fid);
        if numel(siz) == 3
            dose = reshape(dose,siz);
        end
        optResults.dose = dose;

        fid = fopen(weightFileNames{k},'rb');
        weights = fread(fid,'float=>float');
        fclose(fid);
        optResults.weights = weights;
    end

    % load DVH information
    if isfield(optSettings,'prescInfo') & isfield(optSettings.prescInfo,'presc') ...
            & isfield(optSettings.prescInfo.presc,'tissue') ...
            & isfield(optSettings.prescInfo.presc.tissue,'ind')
        optResults.presc = optSettings.prescInfo.presc;

        if readAllOutputFiles == 1
            % calculate a max dose vector
            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} ...
                        = single(dvh(optResults.dose{k},tissMask,dvhbins))';
                end
            end
        else
            dmax = 1.1*max(optResults.dose(:));
            % calculate DVHs and EVHs 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 ...
                    = single(dvh(optResults.dose,tissMask,dvhbins))';
                % calculate the EDVH if presc is non-zero
                if sum(optResults.presc.tissue(m).dMinus) > 0
                    % extract dose voxels
                    doseRatio = zeros(optResults.presc.siz);
                    doseRatio(optResults.presc.tissue(m).ind) = ...
                        optResults.dose(optResults.presc.tissue(m).ind)./ ...
                        optResults.presc.tissue(m).dMinus*100;
                    % remove and Infs and NaNs
                    doseRatio(isinf(doseRatio)) = 0;
                    doseRatio(isnan(doseRatio)) = 0;
                    % calculate edvh
                    edvhbins = [0:Ndvhbins-1]*max(doseRatio(:))/Ndvhbins;
                    optResults.presc.tissue(m).edvhbins = edvhbins;
                    optResults.presc.tissue(m).edvh = ...
                        single(dvh(doseRatio,tissMask,edvhbins))';
                end
            end
        end
    else
        fprintf('Unable to load DVH information: not enough prescription fields available.\n\n');
    end
else
    fprintf('Unable to load beamlet weights and dose distributions: not enough fields available.\n\n');
end

% load the objective function
if isfield(optSettings,'optInfo') & isfield(optSettings.optInfo,'outputFolder') ...
        & isfield(optSettings.optInfo,'objFuncFileName')
    fid = fopen([optSettings.optInfo.outputFolder '/' optSettings.optInfo.objFuncFileName],'rb');
    if fid ~= -1
        optResults.objFunc = fread(fid,'float=>float');
        fclose(fid);
    end
else
    fprintf('Unable to load objective function file: not enough fields available.');
end

% sort all the loaded cell arrays based on batchNames
batch_indices = cellfun(@str2num, optResults.batchNames);
[delme perm_vector] = sort(batch_indices);
optResults.batchNames = optResults.batchNames(perm_vector);
optResults.dose = optResults.dose(perm_vector);
optResults.weights = optResults.weights(perm_vector);

% save file to matlab_files folder
savefile = fullfile(optimizerFolder, 'matlab_files', 'optResults.mat');
save(savefile,'optSettings','optResults');