function optResults = loadOptResults(varargin)
% Loads the weights and dose distribution corresponding to a given input
% filename.  The input filename has the following format:
%
% optResults = loadOptResults(inputFileName) loads all of the files from
%        the linlsq optimization associated with inputFileName.  DVHs are
%        calculated for all of the associated dose files.
% optResults = 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

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

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

% Open the optimization input file
fid = fopen(inputFileName,'r');

% read the file into a cell array, line-by-line
fileText = {};
while 1
    tline = fgetl(fid);
    if ~ischar(tline)
        break; 
    end
    fileText{end+1} = tline;
end
fclose(fid);

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

optResults = [];  % optimization results structure

% search for key fields
for k=1:length(fileText)
    if strcmp(deblank(fileText{k}),'Niterations')
        optResults.optParam.Niterations = str2num(fileText{k+1});
    end
    if strcmp(deblank(fileText{k}),'Nperbatch')
        optResults.optParam.Nperbatch = str2num(fileText{k+1});
    end
    if strcmp(deblank(fileText{k}),'prescription_filename')
        optResults.optParam.prescFileName = fileText{k+1};
    end
    if strcmp(deblank(fileText{k}),'initial_beam_weights_filename')
        optResults.optParam.initialBeamWeightsFileName = fileText{k+1};
    end
    if strcmp(deblank(fileText{k}),'beamlet_header_file')
        optResults.optParam.beamletHeaderFileName = fileText{k+1};
    end
    if strcmp(deblank(fileText{k}),'dose_batch_base_name')
        optResults.optParam.doseBatchBaseName = fileText{k+1};
    end
    if strcmp(deblank(fileText{k}),'dose_batch_extension')
        optResults.optParam.doseBatchExt = fileText{k+1};
    end
    if strcmp(deblank(fileText{k}),'weight_batch_base_name')
        optResults.optParam.weightBatchBaseName = fileText{k+1};
    end
    if strcmp(deblank(fileText{k}),'weight_batch_extension')
        optResults.optParam.weightBatchExt = fileText{k+1};
    end
    if strcmp(deblank(fileText{k}),'obj_func_name')
        optResults.optParam.objFuncName = fileText{k+1};
    end
end

% open the prescription file
fid = fopen([inputFolder '/' optResults.optParam.prescFileName],'r');
% read the entire file into a cell array
fileText = {};
while 1
    tline = fgetl(fid);
    if ~ischar(tline)
        break; 
    end
    fileText{end+1} = tline;
end
fclose(fid);

if ~strcmp(deblank(fileText{1}),'Ntissue')
    error('First line of prescription file must be ''Ntissue''');
else
    Ntissue = str2num(fileText{2});
end

k = 5;  % skip up to the 5th line of the file, which is the tissue number

% read in the tissues
for m=1:Ntissue  % read through the lines for the current tissue
    tissNum = str2num(fileText{k});
    k = k + 2;
    optResults.presc.tissue(m).name = fileText{k};
    k = k + 2;
    optResults.presc.tissue(m).alpha = str2num(fileText{k});
    k = k + 2;
    optResults.presc.tissue(m).betaVPlus = str2num(fileText{k});
    k = k + 2;
    optResults.presc.tissue(m).dVPlus = str2num(fileText{k});
    k = k + 2;
    optResults.presc.tissue(m).vPlus = str2num(fileText{k});
    k = k + 2;
    optResults.presc.tissue(m).betaVMinus = str2num(fileText{k});
    k = k + 2;
    optResults.presc.tissue(m).dVMinus = str2num(fileText{k});
    k = k + 2;
    optResults.presc.tissue(m).vMinus = str2num(fileText{k});
    k = k + 2;
    optResults.presc.tissue(m).betaPlus = str2num(fileText{k});
    k = k + 2;
    optResults.presc.tissue(m).dosePlusFileName = fileText{k};
    k = k + 2;
    optResults.presc.tissue(m).betaMinus = str2num(fileText{k});
    k = k + 2;
    optResults.presc.tissue(m).doseMinusFileName = fileText{k};
    
    % load-in the dPlus inforomation
    fid = fopen([inputFolder '/' optResults.presc.tissue(m).dosePlusFileName],'rb');
    siz = fread(fid,3,'int');  % size of the sparse array
    Nind = fread(fid,1,'int');
    ind = fread(fid,Nind,'int=>int');
    dat = fread(fid,Nind,'float=>float');
    fclose(fid);
    optResults.presc.tissue(m).ind = ind;
    optResults.presc.tissue(m).dPlus = dat;
    
    % load-in the dMinus inforomation
    fid = fopen([inputFolder '/' optResults.presc.tissue(m).doseMinusFileName],'rb');
    siz = fread(fid,3,'int')';  % size of the sparse array
    Nind = fread(fid,1,'int');
    ind = fread(fid,Nind,'int=>int');
    dat = fread(fid,Nind,'float=>float');
    fclose(fid);
    optResults.presc.tissue(m).ind = ind;
    optResults.presc.tissue(m).dMinus = dat;
    
    k = k + 3;  % skip up to the next tissueNum value field
end

optResults.presc.siz = siz;  % size of all of the dose grids

% find which dose and weights files are present
Nbatches = ceil(optResults.optParam.Niterations/optResults.optParam.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 optResults.optParam.doseBatchBaseName ...
        num2str((k-1)*optResults.optParam.Nperbatch) '.' ...
        optResults.optParam.doseBatchExt];
    weightFileNames{k} = [inputFolder optResults.optParam.weightBatchBaseName ...
        num2str((k-1)*optResults.optParam.Nperbatch) '.' ...
        optResults.optParam.weightBatchExt];
    
    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

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,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,siz);
    optResults.dose = dose;
    
    fid = fopen(weightFileNames{k},'rb');
    weights = fread(fid,'float');
    fclose(fid);
    optResults.weights = weights;
end

% 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
    fprintf('Finished calculating DVHs for %s\n',optResults.presc.tissue(m).name);   
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 optResults.optParam.objFuncName],'rb');
optResults.objFunc = fread(fid,'float');
fclose(fid);