function [FIGH,AXH,OBJH,LIKH,ESTH,BH] = plotneigh(d,varargin)
% plotneigh  Plot local behaviour of objective function after estimation.
%
% Syntax
% =======
%
%     [FIGH,AXH,OBJH,LIKH,ESTH,BH] = grfun.plotneigh(D,...)
%
% Input arguments
% ================
%
% * `D` [ struct ] - Structure describing the local behaviour of the
% objective function returned by the [`neighbourhood`](model/neighbourhood)
% function.
%
% Output arguments
% =================
%
% * `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
% ========
%
% * `'plotObj='` [ *`true`* | `false` ] - Plot the local behaviour of the overall
% objective function.
%
% * `'plotLik='` [ *`true`* | `false` ] - Plot the local behaviour of the data
% likelihood component.
%
% * `'plotEst='` [ *`true`* | `false` ] - Mark the actual parameter estimate.
%
% * `'plotBounds='` [ *`true`* | `false` ] - Draw the lower and/or upper bounds
% if they fall within the graph range.
%
% * `'subplot='` [ *`'auto'`* | numeric ] - Subplot division of the figure
% when plotting the results.
%
% * `'linkAxes='` [ `true` | *`false`* ] - Make the vertical axes identical for
% all graphs.
%
% Description
% ============
%
% The data log-likelihood curves are shifted up or down by an arbitrary
% constant to make them fit in the graph; their curvature is preserved.
%
% Example
% ========
%

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

options = passvalopt('grfun.plotneigh',varargin{:});

isplot = options.plotobj || options.plotlik || options.plotest;

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

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

if ~isplot
    return
end

plist = fieldnames(d);
np = numel(plist);

if strcmpi(options.subplot,'auto')
    sub = ceil(sqrt(np));
    if sub*(sub-1) >= np
        sub = [sub-1,sub];
    else
        sub = [sub,sub];
    end
else
    sub = options.subplot;
end

n = prod(sub);
ylim = nan(np,2);

% Force a new figure window.
count = n + 1;

for i = 1 : np
    count = count + 1;
    if count > n
        count = 1;
        FIGH(end+1) = figure(); %#ok<AGROW>
    end
    AXH(end+1) = subplot(sub(1),sub(2),count); %#ok<AGROW>
    x = d.(plist{i}){1};
    xmin = min(x);
    xmax = max(x);
    % Objective func (minus posterior density).
    y1 = d.(plist{i}){2}(:,1);
    % Objective func at optimum.
    x3 = d.(plist{i}){3}(1);
    y3 = d.(plist{i}){3}(2);
    hold('all');
    % Minus log lik of data.
    y2 = d.(plist{i}){2}(:,2);
    [ans,index] = min(abs(x - x3));
    z = y1(index) - y2(index);
    y2 = y2 + z;
    if options.plotobj
        OBJH(end+1) = plot(x,y1); %#ok<AGROW>
    end
    if options.plotlik
        LIKH(end+1) = plot(x,y2); %#ok<AGROW>
    end
    if options.plotest
        ESTH(end+1) = plot(x3,y3, ...
            'marker','*','lineStyle','none','color','red'); %#ok<AGROW>
    end
    if options.plotbounds
        lo = d.(plist{i}){3}(3);
        hi = d.(plist{i}){3}(4);
        if lo >= xmin && lo <= xmax
            BH(end+1) = grfun.vline(lo,'color','red'); %#ok<AGROW>
        end
        if hi >= xmin && hi <= xmax
            BH(end+1) = grfun.vline(hi,'color','red'); %#ok<AGROW>
        end
    end
    grid('on');
    title(plist{i},'interpreter','none');
    axis('tight');
    ylim(i,:) = get(AXH(end),'yLim');
    try %#ok<TRYNC>
        set(AXH(end),'xLim',[xmin,xmax]);
    end
end

if options.linkaxes
    % Sort ylims by the total coverage.
    [ans,index] = sort(ylim*[-1;1]); %#ok<*NOANS,*ASGLU>
    ylim = ylim(index,:);
    linkaxes(AXH,'y');
    % Set ylims to the median coverage.
    set(AXH(end),'yLim',ylim(ceil(np/2),:));
end

end