function[THIS,pstar,objstar,HESS,boundhit,S] = myestimate(THIS,S)
% myestimate  [Not a public function] Estimate parameters.
%
% Backend IRIS function.
% No help provided.

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

% Use fast solve throughout `estimate`. Use regular solve at the end
% with the final parameter estimates to get full solution including
% expansion matrices.
S.fastsolve = true;

% Set Optimization Toolbox options structure.
dooptimoptions();

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

if ~isempty(S.p0)
    
    np = length(S.p0);
    HESS = {zeros(np),zeros(np)};
    % The `assign` and `stdcorr` vectors are used in `objfunc` to reset
    % the model's parameterisation. This is to make sure that the `obj`
    % handle returned as an output of `estimate` will be not be
    % affected by re-scaling the std devs in the output model object.
    % Make sure the model is solved in the very first run.
    if ischar(S.solver)
        if strncmpi(S.solver,'fmin',4)
            if all(isinf(S.pl)) && all(isinf(S.pu))
                [pstar,objstar,EXITFLAG,OUTPUT,grad,HESS{1}] = ...
                    fminunc(@objfunc,S.p0,S.optimset,THIS,S); %#ok<ASGLU>
                LAMBDA = struct('lower',zeros(np,1),'upper',zeros(np,1));
            else
                [pstar,objstar,EXITFLAG,OUTPUT,LAMBDA,grad,HESS{1}] = ...
                    fmincon(@objfunc,S.p0, ...
                    [],[],[],[],S.pl,S.pu,[],S.optimset,THIS,S); %#ok<ASGLU>
            end
        elseif strcmpi(S.solver,'lsqnonlin')
            [pstar,objstar,~,~,~,LAMBDA] = ...
                lsqnonlin(@objfunc,S.p0,S.pl,S.pu,S.optimset,THIS,S);
        end
        % Find lower or upper bound hits.
        boundhit = double(LAMBDA.lower) ~= 0 | double(LAMBDA.upper) ~= 0;
        boundhit = boundhit(:).';
    else
        % User-supplied minimisation routine.
        if isa(S.solver,'function_handle')
            % User supplied function handle.
            f = S.solver;
            args = {};
        else
            % User supplied cell `{func,arg1,arg2,...}`.
            f = S.solver{1};
            args = S.solver(2:end);
        end
        [pstar,objstar,HESS{1}] = ...
            f(@(p) objfunc(p,THIS,S),S.p0,S.pl,S.pu,S.optimset,args{:});
        boundhit = pstar == S.pl | pstar == S.pu;
        boundhit = boundhit(:).';
    end
    
    % Fix numerical inaccuracies when `fmincon` sometimes returns values
    % numerically below lower bounds or above upper bounds.
    dochkbounds();
    
    % Compute contributions of prior distributions to gradient and
    % hessian.
    if any(S.priorindex)
        hessprior = estimateobj.diffprior(S.prior,pstar);
        HESS{2}(:,:) = HESS{2}(:,:) + hessprior;
    end
    
    % Assign estimated parameters, refresh dynamic links, and
    % re-compute steady state and solution. Use regular, not fast,
    % solve to get full solution including also the expansion matrices.
    THIS.Assign = S.assign;
    THIS.stdcorr = S.stdcorr;
    S.fastsolve = false;
    THIS = myupdatemodel(THIS,pstar,S.assignpos,S.stdcorrpos,S);
    
else
    
    % No parameters to be estimated.
    utils.warning('model','No parameters to be estimated.');
    pstar = zeros(1,0);
    objstar = NaN;
    HESS = {zeros(0),zeros(0)};
    boundhit = zeros(1,0);
    
end

% Nested functions.

%**************************************************************************
    function dochkbounds()
        below = pstar < S.pl;
        above = pstar > S.pu;
        if any(below)
            belows = {};
            for ii = find(below)
                belows = [belows,{pstar(ii),S.pl(ii)-pstar(ii),S.plist{ii}}]; %#ok<AGROW>
            end
            utils.warning('model', ...
                ['Final estimate (%g) for this parameter is ', ...
                'numerically below its lower bound by a margin of %g ', ...
                'and will be reset: ''%s''.'], ...
                belows{:});
        end
        if any(above)
            aboves = {};
            for ii = find(above)
                aboves = [aboves,{pstar(ii),pstar(ii)-S.pu(ii),S.plist{ii}}]; %#ok<AGROW>
            end
            utils.warning('model', ...
                ['Final estimate (%g) for this parameter is ', ...
                'numerically above its upper bound by a margin of %g ', ...
                'and will be reset: ''%s''.'], ...
                aboves{:});
        end
        % Reset the out-of-bounds values.
        pstar(below) = S.pl(below);        
        pstar(above) = S.pu(above);
    end
% dochkbounds().

%**************************************************************************
    function dooptimoptions()
        solvername = '';
        if ischar(S.solver)
            solvername = S.solver;
        elseif isa(S.solver,'function_handle')
            solvername = char(S.solver);
        elseif iscell(S.solver)
            solvername = char(S.solver{1});
        end
        switch lower(solvername)
            case 'lsqnonlin'
                algorithm = 'levenberg-marquardt';
                % algorithm = 'trust-region-reflective';
            otherwise
                algorithm = 'active-set';
        end
        oo = {...
            'algorithm',algorithm, ...
            'display',S.display, ...
            'maxIter',S.maxiter, ...
            'maxFunEvals',S.maxfunevals, ...
            'GradObj','off', ...
            'Hessian','off', ...
            'LargeScale','off', ...
            'tolFun',S.tolfun, ...
            'tolX',S.tolx, ...
            };
        if iscell(S.optimset)
            oo = [oo,S.optimset];
        end
        oo(1:2:end) = strrep(oo(1:2:end),'=','');
        S.optimset = optimset(oo{:});
    end
% dooptimoptions().

end