%
%  read a Dakota .out file and parse it.
%
%  [method,dresp,scm,pcm,srcm,prcm]=dakota_out_parse(filei)
%
function [method,dresp,scm,pcm,srcm,prcm]=dakota_out_parse(filei)

if ~nargin
    help dakota_out_parse
    return
end

if ~exist('filei' ,'var') || isempty(filei)
    filei=input('Input file?  ','s');
end
fidi=fopen(sprintf('%s',filei),'r');
if (fidi < 0)
    error('%s could not be opened.',filei);
end

%%  check the first line for the Dakota tabular output file

method=[];
fline=fgetl(fidi);
if ~ischar(fline)
    ieof=1;
    error('File %s is empty.',filei);
end

if strncmpi(fline,'%eval_id',8)
    method='unknown';
    [dresp]=dak_tab_out(fidi,fline);
    return
else
    fseek(fidi,0,'bof');
end

%%  loop through the file to find the Dakota method name

[fline]=findline(fidi,'methodName = ');
if ~ischar(fline)
    return
end

[ntokens,tokens]=fltokens(fline);
method=tokens{1}{3};
display(sprintf('Dakota method=%s.',method));

scm =[];
pcm =[];
srcm=[];
prcm=[];

%%  switch according to the Dakota method

switch lower(method)
%     case {'dot_bfgs','dot_frcg','dot_mmfd','dot_slp','dot_sqp'}
%     case {'npsol_sqp'}
    case {'conmin_frcg','conmin_mfd'}
        [dresp]=opt_any(fidi);
%     case {'optpp_cg','optpp_q_newton','optpp_fd_newton',...
%             'optpp_newton','optpp_pds'}
%     case {'asynch_pattern_search'}
%     case {'coliny_cobyla','coliny_direct','coliny_ea',...
%             'coliny_pattern_search','coliny_solis_wets'}
%     case {'ncsu_direct'}
%     case {'moga','soga'}
%     case {'nl2sol','nlssol_sqp','optpp_g_newton'}
    case {'nond_sampling'}
        [dresp,scm,pcm,srcm,prcm]=nond_samp(fidi);
    case {'nond_local_reliability'}
        [dresp]=nond_locrel(fidi);
%     case {'dace','fsu_quasi_mc','fsu_cvt'}
%     case {'vector_parameter_study','list_parameter_study',...
%             'centered parameter_study','multidim_parameter_study'}
    otherwise
        error('Unrecognized  method ''%s'' in file %s.',...
            method,filei);
end

%%  loop through the file to verify the end

% [fline]=findline(fidi,'<<<<< Single Method Strategy completed');
% if ~ischar(fline)
%     return
% end
display('End of file successfully reached.');
fclose(fidi);

end

%%  function to parse the dakota tabular output file

function [dresp]=dak_tab_out(fidi,fline)

display('Reading Dakota tabular output file.');

%  process column headings of matrix (skipping eval_id)

[ntokens,tokens]=fltokens(fline);
desc=cell (1,ntokens-1);
data=zeros(1,ntokens-1);

for i=1:ntokens-1
    desc(1,i)=cellstr(tokens{1}{i+1});
end
display(sprintf('Number of columns (Dakota V+R)=%d.',ntokens-1));
    
%  process rows of matrix

nrow=0;
while 1
    fline=fgetl(fidi);
    if ~ischar(fline) || isempty(fline)
        break;
    end
    [ntokens,tokens]=fltokens(fline);

%  add row values to matrix (skipping eval_id)

    nrow=nrow+1;
    for i=1:ntokens-1
        data(nrow,i)=tokens{1}{i+1};
    end
end
display(sprintf('Number of rows (Dakota func evals)=%d.',nrow));

%  calculate statistics

%dmean  =mean   (data);
%dstddev=std    (data,0);
[dmean,dstddev,dmeanci,dstddevci]=...
    normfit(data,0.05);

dmin   =min    (data);
dquart1=prctile(data,25);
dmedian=median (data);
dquart3=prctile(data,75);
dmax   =max    (data);

%  same as Dakota scm, Excel correl
dcorrel=corrcoef(data);

%  divide the data into structures for consistency

for i=1:length(desc)
    dresp(i).descriptor=char(desc(i));
    dresp(i).sample    =data(:,i);
    dresp(i).mean      =dmean(i);
    dresp(i).stddev    =dstddev(i);
    dresp(i).meanci    =dmeanci(:,i);
    dresp(i).stddevci  =dstddevci(:,i);
    dresp(i).min       =dmin(i);
    dresp(i).quart1    =dquart1(i);
    dresp(i).median    =dmedian(i);
    dresp(i).quart3    =dquart3(i);
    dresp(i).max       =dmax(i);
end

%  draw box plot

figure
subplot(2,1,1)
plot_boxplot(dresp);

%  draw normal probability plot

subplot(2,1,2)
plot_normplot(dresp);

end

%%  function to parse the nond_sampling output file

function [dresp,scm,pcm,srcm,prcm]=nond_samp(fidi)

ieof=0;

%%  loop through the file to find the Dakota statistics

[fline]=findline(fidi,'Statistics based on ');
if ~ischar(fline)
    return
end

[ntokens,tokens]=fltokens(fline);
nsamp=tokens{1}{4};
display(sprintf('Dakota samples=%d.',nsamp));

%%  loop through the file to find the moments

[fline]=findline(fidi,'Moments for each response function');
if ~ischar(fline)
    return
end

display('Reading moments for each response function:');
ndresp=0;

while ~ieof
    fline=fgetl(fidi);
    if ~ischar(fline)
        ieof=1;
        break;
    end
    if isempty(fline)
        break;
    end
    
%  add new response function and moments

    ndresp=ndresp+1;
    [ntokens,tokens]=fltokens(fline);
    dresp(ndresp).descriptor=tokens{1}{ 1};
    display(sprintf('    %s',dresp(ndresp).descriptor));
    dresp(ndresp).mean      =tokens{1}{ 4};
    dresp(ndresp).stddev    =tokens{1}{ 8};
    dresp(ndresp).coefvar   =tokens{1}{13};
end

display(sprintf('Number of Dakota response functions=%d.',ndresp));

%%  loop through the file to find the confidence intervals

[fline]=findline(fidi,...
    '95% confidence intervals for each response function');
if ~ischar(fline)
    return
end

display('Reading 95% confidence intervals for each response function.');

while ~ieof
    fline=fgetl(fidi);
    if ~ischar(fline)
        ieof=1;
        break;
    end
    if isempty(fline)
        break;
    end
    
%  add confidence intervals to response functions

    [ntokens,tokens]=fltokens(fline);
    for i=1:ndresp
        if strcmpi(tokens{1}{ 1},dresp(i).descriptor)
            dresp(i).meanci  (1,1)=tokens{1}{ 5};
            dresp(i).meanci  (2,1)=tokens{1}{ 6};
            dresp(i).stddevci(1,1)=tokens{1}{12};
            dresp(i).stddevci(2,1)=tokens{1}{13};
            break;
        end
    end
end

%%  loop through the file to find the probabilities

[fline]=findline(fidi,'Probabilities for each response function');
if ~ischar(fline)
    return
end

display('Reading CDF''s for each response function.');

while ~ieof && ~isempty(fline)
    fline=fgetl(fidi);
    if ~ischar(fline)
        ieof=1;
        break;
    end

%  process header line of cdf

    while ~isempty(fline)
        [ntokens,tokens]=fltokens(fline);

%  find response function associated with cdf

        idresp=0;
        for i=1:ndresp
            if strcmpi(tokens{1}{ 6},dresp(i).descriptor)
                idresp=i;
                break;
            end
        end
    
%  skip column headings of cdf

        fline=fgetl(fidi);
        fline=fgetl(fidi);

%  read and add cdf table to response function

        fline=fgetl(fidi);
        icdf=0;
        while ~isempty(fline) && ...
                ~strncmpi(fline,'Cumulative Distribution Function',32)
            [ntokens,tokens]=fltokens(fline);
            icdf=icdf+1;
            dresp(idresp).cdf(icdf,1)=tokens{1}{ 1};
            dresp(idresp).cdf(icdf,2)=tokens{1}{ 2};
            if (ntokens == 4)
                dresp(idresp).cdf(icdf,3)=tokens{1}{ 3};
                dresp(idresp).cdf(icdf,4)=tokens{1}{ 4};
            else
                dresp(idresp).cdf(icdf,3)=NaN;
                dresp(idresp).cdf(icdf,4)=NaN;
            end
            fline=fgetl(fidi);
        end
    end
end

%%  loop through the file to find the scm

[fline]=findline(fidi,'Simple Correlation Matrix');
if ~ischar(fline)
    scm=[];
    return
end

display('Reading simple correlation matrix.');
scm=corrmat(fidi,fline);

%%  loop through the file to find the pcm

[fline]=findline(fidi,'Partial Correlation Matrix');
if ~ischar(fline)
    pcm=[];
    return
end

display('Reading partial correlation matrix.');
pcm=corrmat(fidi,fline);

%%  loop through the file to find the srcm

[fline]=findline(fidi,'Simple Rank Correlation Matrix');
if ~ischar(fline)
    srcm=[];
    return
end

display('Reading simple rank correlation matrix.');
srcm=corrmat(fidi,fline);

%%  loop through the file to find the prcm

[fline]=findline(fidi,'Partial Rank Correlation Matrix');
if ~ischar(fline)
    prcm=[];
    return
end

display('Reading partial rank correlation matrix.');
prcm=corrmat(fidi,fline);

end

%%  function to read a correlation matrix at current file position

function [cmat]=corrmat(fidi,fline)

ieof=0;
cmat.title=fline;

while ~ieof && ~isempty(fline)
    fline=fgetl(fidi);
    if ~ischar(fline)
        ieof=1;
        break;
    end

%  process column headings of matrix

    [ntokens,tokens]=fltokens(fline);
    cmat.column=cell(1,ntokens);
    cmat.row   =cell(1,1);
    cmat.matrix=zeros(1,ntokens);
    
    for i=1:ntokens
        cmat.column(1,i)=cellstr(tokens{1}{i});
    end
    
%  process rows of matrix

    nrow=0;
    while ~ieof
        fline=fgetl(fidi);
        if isempty(fline)
            break;
        end
        [ntokens,tokens]=fltokens(fline);

%  add row heading to matrix

        nrow=nrow+1;
        cmat.row   (nrow,1)=cellstr(tokens{1}{1});

%  add row values to matrix

        for i=2:ntokens
            cmat.matrix(nrow,i-1)=tokens{1}{i};
        end
    end
end

end

%%  function to parse the nond_local_reliability output file

function [dresp]=nond_locrel(fidi)

ieof=0;

%%  loop through the file to find the Dakota statistics

[fline]=findline(fidi,'<<<<< Function evaluation summary');
if ~ischar(fline)
    return
end

[ntokens,tokens]=fltokens(fline);
nfunc=tokens{1}{5};
display(sprintf('Dakota function evaluations=%d.',nfunc));

%%  loop through the file to find the statistics

[fline]=findline(fidi,'MV Statistics for ');
if ~ischar(fline)
    return
end

display('Reading MV statistics for each response function:');
ndresp=0;

while ~ieof && ~isempty(fline) && ...
        strncmpi(fline,'MV Statistics for ',18)

%  add new response function and moments

    ndresp=ndresp+1;
    [ntokens,tokens]=fltokens(fline);
    dresp(ndresp).descriptor=tokens{1}{4};
    display(sprintf('    %s',dresp(ndresp).descriptor));
    fline=fgetl(fidi);
    [ntokens,tokens]=fltokens(fline);
    dresp(ndresp).mean      =tokens{1}{5};
    fline=fgetl(fidi);
    [ntokens,tokens]=fltokens(fline);
    dresp(ndresp).stddev    =tokens{1}{7};

%  read and add importance factors to response function

	idvar=0;
    fline=fgetl(fidi);
    if ~ischar(fline)
        ieof=1;
        break;
    end

    while ~ieof && ~isempty(fline) && ...
            strncmpi(fline,'  Importance Factor for variable ',33)
        [ntokens,tokens]=fltokens(fline);
        idvar=idvar+1;
        dresp(ndresp).desvar(idvar,1)=cellstr(tokens{1}{5});
        dresp(ndresp).impfac(idvar,1)=        tokens{1}{7};

        fline=fgetl(fidi);
        if ~ischar(fline)
            ieof=1;
            break;
        end
    end

%  if importance factors missing, skip to cdf

    if ~idvar
        display('      Importance Factors not available.');
        dresp(ndresp).desvar={};
        dresp(ndresp).impfac=[];
        while ~ieof && ...
                ~strncmpi(fline,'Cumulative Distribution Function',32) && ...
                ~strncmpi(fline,'MV Statistics for ',18) && ...
                ~strncmp (fline,'-',1)
            fline=fgetl(fidi);
            if ~ischar(fline)
                ieof=1;
                break;
            end
        end
    end

%  process header line of cdf

    icdf=0;

    while ~ieof && ~isempty(fline) && ...
            strncmpi(fline,'Cumulative Distribution Function',32)
        [ntokens,tokens]=fltokens(fline);

%  find response function associated with cdf

        idresp=0;
        for i=1:ndresp
            if strcmpi(tokens{1}{6},dresp(i).descriptor)
                idresp=i;
                break;
            end
        end
    
%  skip column headings of cdf

        fline=fgetl(fidi);
        fline=fgetl(fidi);

%  read and add cdf table to response function

        fline=fgetl(fidi);
        while ~isempty(fline) && ...
                ~strncmpi(fline,'MV Statistics for ',18) && ...
                ~strncmp (fline,'-',1)
            [ntokens,tokens]=fltokens(fline);
            icdf=icdf+1;
            dresp(idresp).cdf(icdf,1)=tokens{1}{1};
            dresp(idresp).cdf(icdf,2)=tokens{1}{2};
            if (ntokens == 4)
                dresp(idresp).cdf(icdf,3)=tokens{1}{3};
                dresp(idresp).cdf(icdf,4)=tokens{1}{4};
            else
                dresp(idresp).cdf(icdf,3)=NaN;
                dresp(idresp).cdf(icdf,4)=NaN;
            end
            fline=fgetl(fidi);
        end
    end

%  if cdf missing, skip to end of response function

    if ~icdf
        display('      Cumulative Distribution Function not available.');
        dresp(ndresp).cdf=[];
        while ~ieof && ...
                ~strncmpi(fline,'MV Statistics for ',18) && ...
                ~strncmp (fline,'-',1)
            fline=fgetl(fidi);
            if ~ischar(fline)
                ieof=1;
                break;
            end
        end
    end

end

display(sprintf('Number of Dakota response functions=%d.',ndresp));

end

%%  function to parse any optimization output file

function [dresp]=opt_any(fidi)

%%  search through the file to find the Dakota output data

[fline]=findline(fidi,'<<<<< Iterator');
if ~ischar(fline)
    return
end
display(['  ' deblank(fline)]);

[fline]=findline(fidi,'<<<<< Function evaluation summary');
if ~ischar(fline)
    return
end
display(['  ' deblank(fline)]);

[ntokens,tokens]=fltokens(fline);
nfunc=tokens{1}{5};
% display(sprintf('Dakota function evaluations=%d.',nfunc));

%%  loop through the file to find the best evaluation

[fline]=findline(fidi,'<<<<< Best ');
if ~ischar(fline)
    return
end
dresp.best=[];

display('Reading values for best function evaluation:');

while ischar(fline) && ~isempty(fline) && ...
        strncmpi(fline,'<<<<< Best ',11)
    [ntokens,tokens]=fltokens(fline);

%  read and add best parameter(s)

    if     strncmpi(cellstr(tokens{1}{3}),'parameter', 9)
        display(['  ' deblank(fline)]);

        fline=fgetl(fidi);
	    dresp.best.param     =[];
        dresp.best.descriptor={};

        while ischar(fline) && ~isempty(fline) && ...
                ~strncmpi(fline,'<<<<< Best ',11)
            [ntokens,tokens]=fltokens(fline);
            dresp.best.param     (end+1,1)=        tokens{1}{1};
            dresp.best.descriptor(end+1,1)=cellstr(tokens{1}{2});
            fline=fgetl(fidi);
        end

%  read and add best objective function(s)

    elseif strncmpi(cellstr(tokens{1}{3}),'objective', 9) && ...
           strncmpi(cellstr(tokens{1}{4}),'function' , 8)
        display(['  ' deblank(fline)]);

        fline=fgetl(fidi);
	    dresp.best.of=[];

        while ischar(fline) && ~isempty(fline) && ...
                ~strncmpi(fline,'<<<<< Best ',11)
            [ntokens,tokens]=fltokens(fline);
            dresp.best.of(end+1,1)=        tokens{1}{1};
            fline=fgetl(fidi);
        end

%  read and add best residual term(s)

    elseif strncmpi(cellstr(tokens{1}{3}),'residual', 8) && ...
           strncmpi(cellstr(tokens{1}{4}),'term'    , 4)
        display(['  ' deblank(fline)]);

        fline=fgetl(fidi);
	    dresp.best.res=[];

        while ischar(fline) && ~isempty(fline) && ...
                ~strncmpi(fline,'<<<<< Best ',11)
            [ntokens,tokens]=fltokens(fline);
            dresp.best.res(end+1,1)=        tokens{1}{1};
            fline=fgetl(fidi);
        end

%  read and add best constraint value(s)

    elseif strncmpi(cellstr(tokens{1}{3}),'constraint',10) && ...
           strncmpi(cellstr(tokens{1}{4}),'value'     , 5)
        display(['  ' deblank(fline)]);

        fline=fgetl(fidi);
	    dresp.best.nc=[];

        while ischar(fline) && ~isempty(fline) && ...
                ~strncmpi(fline,'<<<<< Best ',11)
            [ntokens,tokens]=fltokens(fline);
            dresp.best.nc(end+1,1)=        tokens{1}{1};
            fline=fgetl(fidi);
        end

%  read and add best data captured

    elseif strncmpi(cellstr(tokens{1}{3}),'data'    , 4) && ...
           strncmpi(cellstr(tokens{1}{4}),'captured', 8)
        display(['  ' deblank(fline)]);
        [ntokens,tokens]=fltokens(fline);
        dresp.best.eval=        tokens{1}{8};

        fline=fgetl(fidi);

        while ischar(fline) && ~strncmpi(fline,'<<<<< Best ',11)
            fline=fgetl(fidi);
        end

%  read until next best or end

    else
        display(['  ' deblank(fline) '  (not processed)']);

        fline=fgetl(fidi);

        while ischar(fline) && ~strncmpi(fline,'<<<<< Best ',11)
            fline=fgetl(fidi);
        end
    end
end

end

%%  function to find a file line starting with a specified string

function [fline]=findline(fidi,string)

ipos=ftell(fidi);

while 1
    fline=fgetl(fidi);
    if ~ischar(fline)
        break;
    else
        if (strncmpi(fline,string,length(string)))
            return;
        end
    end
end

%  issue warning and reset file position

warning('findline:str_not_found',...
    'String ''%s'' not found in file.',string);
fseek(fidi,ipos,'bof');

end

%%  function to parse a file line into tokens

function [ntokens,tokens]=fltokens(fline)

if ~ischar(fline)
    ntokens=-1;
    tokens={};
    return;
end
if isempty(fline)
    ntokens=0;
    tokens={};
    return;
end

strings=textscan(fline,'%s','delimiter',' :');
%for i=1:length(strings{1})
%    display(sprintf('i=%d; strings{1}{%d}=%s',i,i,strings{1}{i}))
%end
ntokens=0;
tokens{1}{length(strings)}='';

for i=1:length(strings{1})
    if isempty(strings{1}{i})
        continue
    end
    ntokens=ntokens+1;
    inum=sscanf(strings{1}{i},'%f');
    if isempty(inum)
        tokens{1}{ntokens}=strings{1}{i};
%        display(sprintf('i=%d; tokens{1}{%d}=%s',...
%            i,ntokens,tokens{1}{ntokens}))
    else
        tokens{1}{ntokens}=inum;
%        display(sprintf('i=%d; tokens{1}{%d}=%f',...
%            i,ntokens,tokens{1}{ntokens}))
    end
end

end
