function [d,FIGH,AXH,OBJH,LIKH,ESTH,BH] = neighbourhood(this,POS,pct,varargin)
% neighbourhood  Evaluate the local behaviour of the objective function around the estimated parameter values.
%
% Syntax
% =======
%
%     [D,FIGH,AXH,OBJH,LIKH,ESTH,BH] = neighbourhood(M,PS,NEIGH,...)
%
% Input arguments
% ================
%
% * `M` [ model | bkwmodel ] - Model or bkwmodel object.
%
% * `PS` [ poster ] - Posterior simulator (poster) object returned by
% the [`model/estimate`](model/estimate) function.
%
% * `NEIGH` [ numeric ] - The neighbourhood in which the objective function
% will be evaluated defined as multiples of the parameter estimates.
%
% Output arguments
% =================
%
% * `D` [ struct ] - Struct describing the local behaviour of the objective
% function and the data likelihood (minus log likelihood) within the
% specified range around the optimum for each parameter.
%
% The following output arguments are non-empty only if you choose `'plot='`
% true:
%
% * `FIGH` [ numeric ] - Handles to the figures created.
%
% * `AXH` [ numeric ] - Handles to the axes objects created.
%
% * `OBJH` [ numeric ] - Handles to the objective function curves plotted.
%
% * `LIKH` [ numeric ] - Handles to the data likelihood curves plotted.
%
% * `ESTH` [ numeric ] - Handles to the actual estimate marks plotted.
%
% * `BH` [ numeric ] - Handles to the bounds plotted.
%
% Options
% ========
%
% * `'plot='` [ *`true`* | `false` ] - Call the
% [`grfun.plotneigh`](grfun/plotneigh) function from within `neighbourhood`
% to visualise the results.
%
% * `'neighbourhood='` [ struct | *empty* ] - Struct specifying the
% neighbourhood points for some of the parameters; these points will be
% used instead of those based on `NEIGH`.
%
% Plotting options
% =================
%
% See help on [`grfun.plotneigh`](grfun/plotneigh) for options available
% when you choose `'plot='` true.
%
% Description
% ============
%
% In the output database, `D`, each parameter is a 1-by-3 cell array.
% The first cell is a vector of parameter values at which the local
% behaviour was investigated. The second cell is a matrix with two column
% vectors: the values of the overall minimised objective function (as set
% up in the [`estimate`](model/estimate) function), and the values of the
% data likelihood component. The third cell is a vector of four numbers:
% the parameter estimate, the value of the objective function at optimum,
% the lower bound and the upper bound.
%
% Example
% ========
%

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

% Parse required input arguments.
P = inputParser();
P.addRequired('M',@(x) isa(x,'estimateobj'));
P.addRequired('POS',@(x) isa(x,'poster'));
P.addRequired('PCT',@isnumeric);
P.parse(this,POS,pct);

% Parse options.
[options,varargin] = passvalopt('model.neighbourhood',varargin{:});

% Reset persistent variables in the steady state solver.
clear('model/mysstatenonlin');

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

FIGH = [];
AXH = [];
OBJH = [];
LIKH = [];
ESTH = [];
BH = [];

d = struct();
plist = POS.paramList;
pstar = POS.initParam;
np = numel(plist);
npct = numel(pct);
man = options.neighbourhood;

if options.progress
    msg = sprintf('IRIS %s.neighbourhood progress',class(this));
    progress = progressbar(msg);
end

for i = 1 : np
    x = cell(1,3);
    if isstruct(man) && isfield(man,plist{i})
        n = numel(man.(plist{i}));
        pp = cell(1,n);
        pp(:) = {pstar};
        for j = 1 : n
            pp{j}(i) = man.(plist{i})(j);
            x{1}(end+1,1) = pp{j}(i);
        end
    else
        n = npct;
        pp = cell(1,n);
        pp(:) = {pstar};
        for j = 1 : npct
            pp{j}(i) = pstar(i)*pct(j);
            x{1}(end+1,1) = pp{j}(i);
        end
    end
    x{2} = zeros(n,2);
    [x{2}(:,1),x{2}(:,2)] = eval(POS,pp{:}); %#ok<EVLC>
    x{2} = -x{2};
    x{3} = [POS.initParam(i),-POS.initLogPost, ...
        POS.lowerBounds(i),POS.upperBounds(i)];
    d.(plist{i}) = x;
    if options.progress
        update(progress,i/np);
    end
end

if options.plot
    [FIGH,AXH,OBJH,LIKH,ESTH,BH] = grfun.plotneigh(d,varargin{:});
end

end