function [ARG1,ARG2,AR,SGM,FINALCOV,SAVECOUNT] = simulate(THIS,CALLER,DRAW,OPT)
% simulate  [Not a public function] Posterior simulator engine.
%
% Backend IRIS function.
% No help provided.

% -IRIS Toolbox.
% -Copyright (c) 2007-2012 Jaromir Benes & Troy Matheson.

% This function can be called either by `arwm` to draw a chain from the
% posterior distribution, or by `eval` to evaluate the posterior density at
% a particular point.

THETA = [];
LOGPOST = [];
AR = [];
SGM = [];
FINALCOV = [];

%**************************************************************************

npar = length(THIS.paramList);

% Prepackage options.
s = struct();
s.lowerBounds = THIS.lowerBounds(:);
s.upperBounds = THIS.upperBounds(:);
s.lowerBoundsPos = s.lowerBounds > -Inf;
s.upperBoundsPos = s.upperBounds < Inf;
s.lowerBounds = s.lowerBounds(s.lowerBoundsPos);
s.upperBounds = s.upperBounds(s.upperBoundsPos);
s.chkBounds = any(s.lowerBoundsPos) || any(s.upperBoundsPos);
s.isMinusLogPostFunc = isa(THIS.minusLogPostFunc,'function_handle');
if ~s.isMinusLogPostFunc
    s.priorIndex = ...
        cellfun(@(x) isa(x,'function_handle'),THIS.logPriorFunc);
end

% Reset the IRIS steady state solver. This is needed unless the user
% supplied his own solver.
clear('model/mysstatenonlin');

switch CALLER
    case 'arwm'
        doarwm();
        ARG1 = THETA;
        ARG2 = LOGPOST;
    case 'testpar'
        [ARG1,ARG2] = dotestpar();
    case 'eval'
        doeval();
        ARG1 = L;
        ARG2 = K;
end

% Nested functions.

%**************************************************************************
    function doeval()
        % Evaluate log posterior at specified points.
        if ~iscell(DRAW)
            DRAW = {DRAW};
        end
        ndraw = numel(DRAW);
        L = nan(size(DRAW));
        K = nan(size(DRAW));
        for i = 1 : ndraw
            theta = DRAW{i}(:);
            [L(i),K(i)] = dologpostfunc(THIS,theta,s);
        end
    end
% doeval().

%**************************************************************************
    function [t1,t2] = dotestpar()
        % Profile normal for and parfor loops.
        theta = THIS.initParam(:);
        theta = theta(:,ones(1,DRAW));
        L1 = nan(1,DRAW);
        L2 = nan(1,DRAW);
        tic();
        for i = 1 : DRAW
            L1(i) = dologpostfunc(THIS,theta(:,i),s);
        end
        t1 = toc();
        tic();
        parfor (i = 1 : DRAW,DRAW)
            L2(i) = dologpostfunc(THIS,theta(:,i),s); %#ok<PFOUS>
        end
        t2 = toc();
    end
% doeval().

%**************************************************************************
    function doarwm()
        
        % Adaptive random walk Metropolis simulator.
        nalloc = min(DRAW,OPT.saveevery);
        issave = OPT.saveevery <= DRAW;
        if issave
            dochksaveoptions();
        end
        
        sgm = OPT.initscale;
        if OPT.burnin < 1
            burnin = round(OPT.burnin*DRAW);
        else
            burnin = OPT.burnin;
        end
        
        ndraw = DRAW + burnin;
        gamma = OPT.gamma;
        k1 = OPT.adaptscale;
        k2 = OPT.adaptproposalcov;
        targetAR = OPT.targetar;
        
        isAdaptiveScale = ~isnan(gamma) && ~isinf(gamma) && k1 > 0;
        isAdaptiveShape = ~isnan(gamma) && ~isinf(gamma) && k2 > 0;
        isAdaptive = isAdaptiveScale || isAdaptiveShape;
        
        theta = THIS.initParam(:);
        P = chol(THIS.initProposalCov).';
        logPost = dologpostfunc(THIS,theta,s);
        
        % Pre-allocate output data.
        THETA = zeros(npar,nalloc);
        LOGPOST = zeros(1,nalloc);
        AR = zeros(1,nalloc);
        if isAdaptiveScale
            SGM = zeros(1,nalloc);
        else
            SGM = sgm;
        end
        
        if OPT.progress
            progress = progressbar('IRIS poster.arwm progress');
        elseif OPT.esttime
            eta = esttime('IRIS poster.arwm is running');
        end
        
        % Main loop.
        naccepted = 0;
        count = 0;
        SAVECOUNT = 0;
        for j = 1 : ndraw
            % Propose a new theta.
            u = randn(npar,1);
            newTheta = theta + sgm*P*u;
            newLogPost = dologpostfunc(THIS,newTheta,s);
            % Decide if we accept the new theta.
            alpha = min(1,exp(newLogPost-logPost));
            accepted = rand() <= alpha;
            if accepted
                logPost = newLogPost;
                theta = newTheta;
            end
            % Adapt the scale and/or proposal covariance.
            if isAdaptive
                nu = j^(-gamma);
                phi = nu*(alpha - targetAR);
                if isAdaptiveScale
                    phi1 = k1*phi;
                    sgm = exp(log(sgm) + phi1);
                end
                if isAdaptiveShape
                    phi2 = k2*phi;
                    unorm2 = u.'*u;
                    z = sqrt(phi2/unorm2)*u;
                    P = cholupdate(P.',P*z).';
                end
            end
            % Add the j-th theta to the chain unless it's burn-in sample.
            if j > burnin
                count = count + 1;
                naccepted = naccepted + double(accepted);
                % Paremeter draws.
                THETA(:,count) = theta;
                % Value of log posterior at the current draw.
                LOGPOST(count) = logPost;
                % Acceptance ratio so far.
                AR(count) = naccepted / (j-burnin);
                % Adaptive scale factor.
                if isAdaptiveScale
                    SGM(count) = sgm;
                end
                if count == OPT.saveevery ...
                        || (issave && j == ndraw)
                    count = 0;
                    dosave();
                end
            end
            % Update the progress bar or estimated time.
            if OPT.progress
                update(progress,j/ndraw);
            elseif OPT.esttime
                update(eta,j/ndraw);
            end
        end
        
        FINALCOV = P*P.';
        
        if issave
            % Save master file with the following information
            % * `SAVECOUNT` -- the total number of files;
            % * `DRAW` -- the total number of draws;
            PLIST = THIS.paramList; %#ok<NASGU>
            save(OPT.saveas,'PLIST','SAVECOUNT','DRAW');
        end
        
        function dochksaveoptions()
            if isempty(OPT.saveas)
                utils.error('poster', ...
                    'Proper file name must be specified for SAVEAS.');
            end
            [p,t] = fileparts(OPT.saveas);
            OPT.saveas = fullfile(p,t);
        end
        
        function dosave()
            SAVECOUNT = SAVECOUNT + 1;
            filename = [OPT.saveas,sprintf('%g',SAVECOUNT)];
            save(filename,'THETA','LOGPOST','-v7.3');
            togo = ndraw-j;
            if togo == 0
                THETA = [];
                LOGPOST = [];
                AR = [];
                SGM = [];
            elseif togo < nalloc
                THETA = THETA(:,1:togo);
                LOGPOST = LOGPOST(1:togo);
                AR = AR(1:togo);
                if isAdaptiveScale
                    SGM = SGM(1:togo);
                end
            end
        end
        
    end
% doarwm().

end

%**************************************************************************
function [x,K] = dologpostfunc(this,p,s)
% dologpostfunc  Evalute posterior density for given parameters.
% This is a subfunction, and not a nested function, so that we can later
% implement a parfor loop (parfor does not work with nested functions).

K = NaN;
if ~s.chkBounds || ...
        (all(p(s.lowerBoundsPos) >= s.lowerBounds) ...
        && all(p(s.upperBoundsPos) <= s.upperBounds))
    if s.isMinusLogPostFunc
        % Evaluate log posterior.
        x = -this.minusLogPostFunc(p,this.minusLogPostFuncArgs{:});
    else
        % Sum up log likelihood and priors.
        x = 0;
        % Evaluate priors.
        for k = find(s.priorIndex)
            x = x + this.logPriorFunc{k}(p(k));
            if isinf(x)
                return
            end
        end
        % Evaluate minus log likelihood.
        K = -this.minusLogLikFunc(p,this.minusLogLikFuncArgs{:});
        x = x + K;
    end
else
    % Out of bounds.
    x = -Inf;
end

end
% dologpostfunc().