function varargout = soda(varargin)
%
% <a href="matlab:web(fullfile(sodaroot,'html','soda.html'),'-helpbrowser')">View HTML documentation for this function in the help browser</a>
%
% If you haven't used the SODA documentation before, you need to
% install it by running >> soda -docinstall
% <a href="matlab:soda -docinstall">Install documentation now</a>
%
%
% revision of SODA algorithm by J.A. Vrugt
%
% Author         : Jurriaan H. Spaaks
% Date           : August 2007
% Matlab release : 2007a on Win32


if ischar(varargin{1})
    
    sodaInitialize(varargin)

elseif isstruct(varargin{1})

    % display disclaimer:
    disp_disclaimer

    sodaPar = varargin{1};

    % verify the integrity of sodaPar's fieldnames:
    sodaVerifyFieldNames(sodaPar,'check')

    % % These parameters can be derived from the existing settings:
    % Number of parameters to be optimized:
    sodaPar.nOptPars = numel(sodaPar.parNames);
    sodaPar.nStatesKF = numel(sodaPar.stateNamesKF);
    sodaPar.nStatesNOKF = numel(sodaPar.stateNamesNOKF);
    sodaPar.nObjs = numel(sodaPar.objCallStrs);

    % load default settings from file:
    if isfield(sodaPar,'useIniFile')
        fileStr = fullfile(sodaroot,sodaPar.useIniFile);
        sodaPar = sodaLoadSettings(sodaPar,fileStr);
    else
        fileStr = fullfile(sodaroot,'soda-default-settings.ini');
        sodaPar = sodaLoadSettings(sodaPar,fileStr);
    end
    clear fileStr

    % initialize the uniform random generator:
    try
        rand('twister',sodaPar.randSeed);
    catch
        if sodaPar.verboseOutput
            disp(['Using the rand seed method instead of the ',...
                char(10),'default twister method. ',char(10),...
                '<a href="matlab:doc rand">doc rand</a>'])
        end
        rand('seed',sodaPar.randSeed);
    end


    % initialize the Gaussian random generator:
    try
        randn('state',sodaPar.randSeed);
    catch
        if sodaPar.verboseOutput
            disp(['Using the randn state method instead of the ',...
                char(10),'default twister method. ',char(10),...
                '<a href="matlab:doc rand">doc rand</a>'])
        end
        randn('seed',sodaPar.randSeed);
    end



    % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
    % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
    % % % % % %                                             % % % % % %
    % % % % % %      SODA  INITIALIZATION FINISHED       % % % % % %
    % % % % % %                                             % % % % % %
    % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
    % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %


    % check if any inconsistencies can be identified from sodaPar's fields:
    sodaPar = check_input_integrity(sodaPar,nargout);

    % number of samples per complex:
    sodaPar.nSamplesPerCompl = sodaPar.nSamples/sodaPar.nCompl;

    %number of offspring per complex:
    sodaPar.nOffspringPerCompl = sodaPar.nOffspring/sodaPar.nCompl;

    % Calculate the parameters in the exponential power density function:
    [sodaPar.cBeta,sodaPar.omegaBeta] = sodaCalcCbWb(sodaPar);
    
    % specify which states are visualized:
    if strncmp(sodaPar.modeStr,'soda',4) && ~isfield(sodaPar,'nStatesViz')
        sodaPar.nStatesViz = 1:sodaPar.nStatesKF;
    end
    
    
    % Indicate the meaning of each column in 'evalResults':
    sodaPar.iterCol = 1;                                     % model evaluation counter
    sodaPar.parCols = sodaPar.iterCol+(1:sodaPar.nOptPars);  % model parameters
    sodaPar.objCols = sodaPar.parCols(end)+(1:sodaPar.nObjs);% objective scores
    sodaPar.paretoCol = sodaPar.objCols(end)+1;              % pareto assimilation of the objective scores 
    sodaPar.evalCol = sodaPar.paretoCol + 1;                 % whether the model has been evaluated with
                                                             % the parameters in the current row
    
    sodaPar.optStartTime = datestr(now,31);
    
    %initialize the gelman-rubin statistic record:
    critGelRub = [];


    sodaStartTime = now;
    disp([char(10),upper(sodaPar.modeStr) ' run started on: ',datestr(sodaStartTime,...
        'mmmm dd, yyyy'),' ',datestr(sodaStartTime,'HH:MM:ss')])

    if sodaPar.startFromUniform

        switch sodaPar.sampleDrawMode
            case 'stratified random'
                % Startified uniform random sampling of points
                % in parameter space:
                randomDraw = sodaStratRandDraw(sodaPar);
            case 'stratified'
                % Startified sampling of points in parameter space:
                randomDraw = sodaStratDraw(sodaPar);
        end

        evalResults = [[1:sodaPar.nSamples]',randomDraw,...
            repmat([repmat(NaN,[1,sodaPar.nObjs+1]),false],[sodaPar.nSamples,1])];

        evalResults = sodaCalcObjScore(sodaPar,evalResults);
        evalResults(:,sodaPar.paretoCol) = -sodaCalcPareto(evalResults(:,sodaPar.objCols),sodaPar.paretoMethod);
        nModelEvals = size(evalResults,1);

        % partition 'evalResults' into complexes and initialize the sequences:
        complexes = sodaPartComplexes(sodaPar,evalResults);
        sequences = complexes;
        metropolisRejects = repmat(NaN,size(sequences));

    else
        % continue with an earlier run
        reply = input([char(10),'Please provide the filename of the *.mat ',...
            'file containing',char(10),'variables ',...
            char(39),'critGelRub',char(39),', ',...
            char(39),'metropolisRejects',char(39),', ',...            
            char(39),'evalResults',char(39),', and ',...
            char(39),'sequences',char(39),'.',char(10),'>> ']);
        if exist(reply,'file')
            eval(['load(',char(39),reply,char(39),',',...
                char(39),'critGelRub',char(39),',',...
                char(39),'evalResults',char(39),',',...
                char(39),'metropolisRejects',char(39),',',...
                char(39),'sequences',char(39),')'])
        else
            error('No such file exists on the path.')
        end

        nRecords = max(evalResults(:,sodaPar.iterCol));
        while true
            reply2 = input([char(10),'Would you like to re-order the rows of ',...
                'variable ',char(39),'evalResults',char(39),'?',char(10),...
                char(39),'y',char(39),'  yes',char(10),...
                char(39),'n',char(39),'  no',char(10),...            
                char(39),'a',char(39),'  abort',char(10),...
                '>> ']);
            if strcmp(reply2,'y')
                disp(['Reordering rows of ',char(39),'evalResults',char(39),'...'])
                TMP = sortrows(evalResults,sodaPar.objCol);
                evalResults = [(1:size(TMP,1))',TMP(:,2:sodaPar.evalCol)];
                clear TMP
                disp([char(8),'Done'])
                break
            elseif strcmp(reply2,'n')
                disp('Not re-ordering.')
                break                
            elseif strcmp(reply2,'a')
                error('User abort.')
            else
                disp('You used an unknown command input.')
            end
        end

        % partition 'evalResults' into complexes:
        complexes = sodaPartComplexes(sodaPar,evalResults);
        sodaPar.converged = sodaAbortOptim(critGelRub,sodaPar);
        % initialize the sequences:
        if isempty(sequences)
            sequences = complexes;
            nModelEvals = max(evalResults(:,sodaPar.iterCol));
        else
            nTrownAway = mod(nRecords-sodaPar.nSamples,sodaPar.nOffspring);
            nModelEvals = nRecords-nTrownAway;
            if sodaPar.verboseOutput
                disp(['Throwing away last ',num2str(nTrownAway), ' model evaluations.'])
                disp([upper(sodaPar.modeStr),' run continues at ',num2str(nModelEvals+1), ' model evaluations.'])
            end
        end

        
        % temporarily reset the value of 'sodaPar.startFromUniform' in
        % order to invoke an additional test in 'check_input_integrity'
        sodaPar.startFromUniform = true;
        check_input_integrity(sodaPar,nargout)
        % restore the old value
        sodaPar.startFromUniform = false;
        
        if sodaPar.doPlot
            eval(sodaPar.visualizationCall)
        end

    end

    % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
    % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
    % % % % % %                                             % % % % % %
    % % % % % %           SODA  START OF MAIN LOOP          % % % % % %
    % % % % % %                                             % % % % % %
    % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
    % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %

    sodaEstimateAcceptance(sodaPar,evalResults);

    try
        while sodaContinue(sodaPar,nModelEvals)

            for iCompl=1:sodaPar.nCompl

                % select the current complex:
                curCompl = complexes(:,:,iCompl);
                
                % select the current sequence:
                curSeq = sequences(:,:,iCompl);
                
                for iOffspring = 1:sodaPar.nOffspringPerCompl

                    % propose offspring based on statistical properties of the
                    % current complex and sequence:
                    nModelEvals = nModelEvals + 1;
                    propChildren(iOffspring,:,iCompl) = sodaGenerateOffspring(sodaPar,curCompl,curSeq);
                    propChildren(iOffspring,sodaPar.iterCol,iCompl) = nModelEvals;
                    
                end % iOffspring = 1:sodaPar.nOffspringPerCompl

            end % iCompl=1:sodaPar.nCompl

            % calculate objective for unevaluated rows of evalResults:
            propChildren = sodaCalcObjScore(sodaPar,propChildren);
            
            % recalculate the global pareto scores based on the objectives
            % scores for all points evaluated so far:
            [sequences,complexes,propChildren] = sodaRecalcPareto(sodaPar,sequences,complexes,propChildren);
            
            sequences = sodaPrepSeqArray(sodaPar,sequences);
            metropolisRejects = sodaPrepRejArray(sodaPar,metropolisRejects);
            
           
            for iCompl=1:sodaPar.nCompl

                % select the current complex:
                curCompl = complexes(:,:,iCompl);

                % select the current sequence:
                curSeq = sequences(:,:,iCompl);
                
                for iOffspring = 1:sodaPar.nOffspringPerCompl

                    % evolve the current complex by accepting the proposed child...
                    % or by re-accepting the last child from the current sequence:
                    propChild = propChildren(iOffspring,:,iCompl);
                    [curCompl,acceptedChild] = sodaEvolveComplex(sodaPar,curCompl,curSeq,propChild);

                    rNumber = size(sequences,1)-sodaPar.nOffspringPerCompl+iOffspring;

                    % add the accepted child to the 'sequences' array:
                    sequences(rNumber,:,iCompl) = acceptedChild;

                    % add the accepted child to the 'evalResults' array:
                    evalResults = [evalResults;acceptedChild];
                    
                    if isequal(propChild,acceptedChild)
                        metropolisRejects(rNumber,:,iCompl) = repmat(NaN,[1,sodaPar.evalCol]);
                    else
                        metropolisRejects(rNumber,:,iCompl) = propChild;
                    end
                end % iOffspring = 1:sodaPar.nOffspringPerCompl

            end % iCompl=1:sodaPar.nCompl

            
            % recalculate pareto scores based on objective scores in
            % evalResults:
            evalResults(:,sodaPar.paretoCol) = -sodaCalcPareto(evalResults(:,sodaPar.objCols));
            
            % partition 'evalResults' into complexes:
            complexes = sodaPartComplexes(sodaPar,evalResults);

            iGeneration = (nModelEvals-sodaPar.nSamples)/sodaPar.nOffspring;

            %determine the Gelman-Rubin scale reduction (convergence) statistic:
            critGelRub(iGeneration,:) = [nModelEvals,sodaGelmanRubin(sodaPar,sequences)];

            % determine whether all parameters meet the Gelman-Rubin scale
            % reduction criterion
            sodaPar.converged = sodaAbortOptim(critGelRub,sodaPar);

            if sodaPar.doPlot &&...
                    (mod(iGeneration,sodaPar.drawInterval)==0||...
                    iGeneration==1)

                eval(sodaPar.visualizationCall)

            end % sodaPar.doPlot &&...

            if mod(iGeneration,sodaPar.saveInterval)==0||...
                    iGeneration==1

                save([sodaPar.modeStr,'-tempstate.mat'])

            end %  mod(iGeneration,...


        end % sodaContinue(sodaPar,nModelEvals)

        % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
        % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
        % % % % % %                                             % % % % % %
        % % % % % %          SODA END OF MAIN LOOP              % % % % % %
        % % % % % %                                             % % % % % %
        % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
        % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %

        sodaPar.optEndTime = now;
        nDaysPerParSet = (datenum(sodaPar.optEndTime)-datenum(sodaPar.optStartTime))/size(evalResults,1);
        if nDaysPerParSet<(1/(24*60))
            % measure in seconds
            v = nDaysPerParSet*24*60*60;
            u = 'seconds';
        elseif nDaysPerParSet<(1/24)
            % measure in minutes
            v = nDaysPerParSet*24*60;
            u = 'minutes';
        elseif nDaysPerParSet<1
            % measure in hours
            v = nDaysPerParSet*24;
            u = 'hours';
        end
        sodaPar.avgWalltimePerParSet = [sprintf('%f',v),' ',u];

        disp([sodaPar.modeStr,' run completed on: ',datestr(sodaPar.optEndTime,...
            'mmmm dd, yyyy'),' ',datestr(sodaPar.optEndTime,'HH:MM:ss')])

        varargout = prepOutput(nargout,evalResults,critGelRub,sequences,metropolisRejects,sodaPar);

    catch
        
        % there is still a minor bug here somewhere, since the software
        % does not seem to ever reach the catch part.

        varargout = prepOutput(nargout,evalResults,critGelRub,sequences,metropolisRejects,sodaPar);

    end

end








% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% % % % % %                                             % % % % % %
% % % % % %       SODA LOCAL FUNCTIONS START HERE       % % % % % %
% % % % % %                                             % % % % % %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %





% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
function sodaPar = check_input_integrity(sodaPar,nOut)


if all(nOut<[2:5])
    error(['Function ',39,mfilename,39,' should have at least 2 output arguments.',char(10),...
        '[evalResults,critGelRub]=soda(sodaPar)',char(10),...
        '[evalResults,critGelRub,sodaPar]=soda(sodaPar)',char(10),...
        '[evalResults,critGelRub,sequences,sodaPar]=soda(sodaPar)',char(10),...
        '[evalResults,critGelRub,sequences,sodaPar,bayesStore]=soda(sodaPar)',char(10)])
elseif all(nOut>[2:5])
    error(['Function ',39,mfilename,39,' should have at most 5 output arguments.',char(10),...
        '[evalResults,critGelRub]=soda(sodaPar)',char(10),...
        '[evalResults,critGelRub,sodaPar]=soda(sodaPar)',char(10),...
        '[evalResults,critGelRub,sequences,sodaPar]=soda(sodaPar)',char(10),...
        '[evalResults,critGelRub,sequences,sodaPar,bayesStore]=soda(sodaPar)',char(10)])
end

if sodaPar.nCompl<2
    error(['In order to determine the Gelman-Rubin convergence',10,...
        'criterion, at least 2 complexes should be used.'])
end

x = sodaPar.nSamples/sodaPar.nCompl;
if round(x)~=x
    error(['Unable to uniformly distribute number of samples (',...
        num2str(sodaPar.nSamples),') over the given number ',10,...
        'of complexes (',num2str(sodaPar.nCompl),').'])
elseif x<=0
    error(['Variable ',39,'sodaPar.nSamples',39,' must be larger than zero.'])
end

x = sodaPar.nOffspring/sodaPar.nCompl;
if round(x)~=x
    error(['Unable to uniformly distribute number of offspring (',...
        num2str(sodaPar.nOffspring),') over the given number ',10,...
        'of complexes (',num2str(sodaPar.nCompl),').'])
elseif x<=0
    error(['Variable ',39,'sodaPar.nOffspring',39,' must be larger than zero.'])
end
clear x

if sodaPar.convUseLastFraction<=0
    error(['Value of parameter ',39,'sodaPar.convUseLastFraction',39,' should be larger than 0.'])
end

if sodaPar.convUseLastFraction>1
    error(['Value of parameter ',39,'sodaPar.convUseLastFraction',39,' should be smaller than or equal to 1.'])
end

if sodaPar.startFromUniform &&...
        mod((sodaPar.nModelEvalsMax-sodaPar.nSamples),sodaPar.nOffspring)~=0
    
    suggestedValue = sodaPar.nSamples+ceil((sodaPar.nModelEvalsMax-...
        sodaPar.nSamples)/sodaPar.nOffspring)*sodaPar.nOffspring;
    
    error(['Generating ',char(39),'sodaPar.nOffspring',char(39),...
        ' (',num2str(sodaPar.nOffspring),') descendants per ',...
        'generation will not ',char(10),'yield exactly ',char(39),...
        'sodaPar.nModelEvalsMax',char(39),' (',...
        num2str(sodaPar.nModelEvalsMax),') model evaluations ',...
        'for any',char(10),'integer number of generations, given ',...
        'that ',char(39),'sodaPar.nSamples',char(39),' equals ',...
        num2str(sodaPar.nSamples),'.',char(10),'Suggested value = ',num2str(suggestedValue),'.'])
    

end

if ~all(sodaPar.parSpaceLoBound<sodaPar.parSpaceHiBound)
    error(['Parameter domain boundaries incorrectly specified. Check parameters ',...
        char(39),'sodaPar.parSpaceLoBound',char([39,10]),...
        'and ',char(39),'sodaPar.parSpaceHiBound',char(39),'.'])
end

switch sodaPar.modeStr
    case 'scemua'
        rmFields = {'errModelCallStr';...
                    'measErr';...
                    'stateSpaceHiBound';...
                    'stateSpaceLoBound';...
                    'stochForcePars';...
                    'stochForceParsIx'};
        
        for k=1:size(rmFields,1)
            if isfield(sodaPar,rmFields{k,1})
                disp(['Removing field ',char(39),rmFields{k,1},char(39),' from sodaPar.'])
                sodaPar = rmfield(sodaPar,rmFields{k,1});
            end
        end
        clear k
        clear rmFields
    case 'soda'
    case 'reset'
        rmFields = {'errModelCallStr';...
                    'stochForcePars';...
                    'stochForceParsIx'};
        
        for k=1:size(rmFields,1)
            if isfield(sodaPar,rmFields{k,1})
                disp(['Removing field ',char(39),rmFields{k,1},char(39),' from sodaPar.'])
                sodaPar = rmfield(sodaPar,rmFields{k,1});
            end
        end
        clear k
        clear rmFields
        
end



if strncmp(sodaPar.modeStr,'soda',4) && ~all(sodaPar.stateSpaceLoBound<sodaPar.stateSpaceHiBound)
    error(['State domain boundaries incorrectly specified. Check parameters ',...
        char(39),'sodaPar.stateSpaceLoBound',char([39,10]),...
        'and ',char(39),'sodaPar.stateSpaceHiBound',char(39),'.'])
end


if strcmp(sodaPar.optMethod,'error minimization') &...
        ~isfield(sodaPar,'nMeasurements')
    error(['Variable ',char(39),'sodaPar',char(39),...
        ' must have field ',char(39),'nMeasurements',...
        char(39),' when',char(10),'used in ',char(39),'error minimization',...
        char(39),' mode.'])
end

if ~strcmp(sodaPar.optMethod,{'error minimization','likelihood','direct probability'})
    error(['Unknown optimization type in ',char(39),'sodaPar.optMethod',char(39),'.'])
end


switch sodaPar.modeStr
    case 'scemua'
        shouldBeRowList = {'parNames','parSpaceLoBound','parSpaceHiBound'};
    case {'soda','reset'}
        shouldBeRowList = {'parNames','parSpaceLoBound','parSpaceHiBound','stateNamesKF',...
        'stateNamesNOKF','stateSpaceLoBound','stateSpaceHiBound'};
end
for k=1:numel(shouldBeRowList)
    if isfield(sodaPar,shouldBeRowList{k})
        eval(['TMP = sodaPar.',shouldBeRowList{k},';'])
        if size(TMP,1)>1
            error('SODA:shouldBeRowVariable',['Option ',char(39),'sodaPar.',...
                shouldBeRowList{k},char(39),' should be a 1xN variable.'])
        end
    end
end
clear shouldBeRowList TMP k



if strncmp(sodaPar.modeStr,'soda',4) && ~all(ismember(sodaPar.obsPoints{1,2},sodaPar.simPoints{1,2}))
   error(['All elements of ',char(39),'sodaPar.obsPoints',char(39),...
       ' should occur in ',char(39),'sodaPar.simPoints',char(39),'.'])
end

C ={'doPlot','realPartOnly','startFromUniform',...
    'verboseOutput','plotEnsemble'};
for k=1:numel(C)
    if isfield(sodaPar,C{k})
        str = ['IO=~islogical(sodaPar.',C{k},');'];
        eval(str)
        if IO
            error(['Field ',char(39),C{k},char(39), ' should be logical.'])
        end
    end
end
clear C



switch sodaPar.modeStr
    case {'scemua'}
        if sodaPar.nEnsembleMembers~=1
            disp(['Resetting value of ',char(39),'sodaPar.nEnsembleMembers',char(39),' to 1.'])
            sodaPar.nEnsembleMembers = 1;
        end
        if sodaPar.forceRecalcFraction~=0
            disp(['Resetting value of ',char(39),'sodaPar.forceRecalcFraction',char(39),' to 0.'])
            sodaPar.forceRecalcFraction = 0;
        end
    case 'soda'
        % nothing to do here.
    case {'reset'}
        if sodaPar.nEnsembleMembers~=1
            disp(['Resetting value of ',char(39),'sodaPar.nEnsembleMembers',char(39),' to 1.'])
            sodaPar.nEnsembleMembers = 1;
        end
        if sodaPar.forceRecalcFraction~=0
            disp(['Resetting value of ',char(39),'sodaPar.forceRecalcFraction',char(39),' to 0.'])
            sodaPar.forceRecalcFraction = 0;
        end
        if sodaPar.measErr~=0
            disp(['Resetting value of ',char(39),'sodaPar.measErr',char(39),' to 0.'])
            sodaPar.measErr = 0;
        end
        
end




% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
function disp_disclaimer

disp([10,...
    '% This is a pre-alpha release of the SODA simultaneous ',10,...
    '% parameter optimization and data assimilation software.',10,10,...
    '% See <a href=',34,'matlab:web(fullfile(sodaroot,',...
    39,'html',39,',',39,'gpl.txt',39,'),',39,...
    '-helpbrowser',39,')',34,'>link</a> for license.',10,10])




function outCell = prepOutput(nOut,evalResults,critGelRub,sequences,metropolisRejects,sodaPar)

evalResults = sortrows(evalResults,sodaPar.iterCol);

if nOut==2
    outCell{1}=evalResults;
    outCell{2}=critGelRub;
elseif nOut==3
    outCell{1}=evalResults;
    outCell{2}=critGelRub;
    outCell{3}=orderfields(sodaPar);
elseif nOut==4
    outCell{1}=evalResults;
    outCell{2}=critGelRub;
    outCell{3}=sequences;
    outCell{4}=orderfields(sodaPar);
elseif nOut==5
    outCell{1}=evalResults;
    outCell{2}=critGelRub;
    outCell{3}=sequences;
    outCell{4}=metropolisRejects;
    outCell{5}=orderfields(sodaPar);
else
    error(['Function ',39,mfilename,39,' should have at least 2 output arguments.'])
end
