%
%  plot a relative histogram along with the normal distribution.
%
%  []=plot_hist_norm(rfunc1   ,rfunc2  ,hmin,hmax,hnint)
%  []=plot_hist_norm(data,desc,mu,sigma,hmin,hmax,hnint)
%
function []=plot_hist_norm(varargin)

if ~nargin
    help plot_hist_norm
    return
end

%%  process input data and assemble into matrices as needed

ivar=1;
if isstruct(varargin{ivar})
    rfunc1=varargin{ivar};
    ivar=ivar+1;
    
    desc=cell (1,length(rfunc1));
    ldata=zeros(1,length(rfunc1));
    for i=1:length(rfunc1)
        ldata(i)=length(rfunc1(i).sample);
    end
    data=zeros(max(ldata),length(rfunc1));
    data(:,:)=NaN;

    for i=1:length(rfunc1)
        desc(i)=cellstr(rfunc1(i).descriptor);
        data(1:ldata(i),i)=rfunc1(i).sample;
    end
else
    data=varargin{ivar};
    ivar=ivar+1;
    
    ldata(1:size(data,2))=size(data,1);

    if ivar <= nargin && iscell(varargin{ivar})
        desc=varargin{ivar};
        ivar=ivar+1;
    else
        desc={};
    end 
end

if     ivar <= nargin && isstruct(varargin{ivar})
    rfunc2=varargin{ivar};
    ivar=ivar+1;

    mu   =zeros(1,length(rfunc2));
    sigma=zeros(1,length(rfunc2));

    for i=1:length(rfunc2)
        mu   (i)=rfunc2(i).mean;
        sigma(i)=rfunc2(i).stddev;
    end
elseif ivar+1 <= nargin
    if isnumeric(varargin{ivar})
        if ~isempty(varargin{ivar})
            mu   =varargin{ivar};
        else
            mu   =mean(data);
            display('Using calculated means.')
        end
    end
    ivar=ivar+1;
    if isnumeric(varargin{ivar})
        if ~isempty(varargin{ivar})
            sigma=varargin{ivar};
        else
            sigma=std(data);
            display('Using calculated standard deviations.')
        end
    end
    ivar=ivar+1;
end

%%  generate the intervals

if ivar <= nargin && ~isempty(varargin{ivar})
    hmin=varargin{ivar};
else
    hmin=min(min(data));
end
ivar=ivar+1;
if ivar <= nargin && ~isempty(varargin{ivar})
    hmax=varargin{ivar};
else
    hmax=max(max(data));
end
ivar=ivar+1;
if ivar <= nargin && ~isempty(varargin{ivar})
    hnint=varargin{ivar};
else
    hnint=50;
end
ivar=ivar+1;
edges=hmin:(hmax-hmin)/hnint:hmax;

%%  generate the histogram counts and make them relative

%  note that for the histc function:
%  n(k) counts the value x(i) if edges(k) <= x(i) < edges(k+1).
%  The last bin counts any values of x that match edges(end).
%  Values outside the values in edges are not counted.
%  Use -inf and inf in edges to include all non-NaN values.

dhistc=histc(data,edges);
for i=1:size(data,2)
    dbelow(i)  =length(find(data(:,i)<edges(  1)))/ldata(i);
    dhistc(:,i)=dhistc(:,i)                       /ldata(i);
    dabove(i)  =length(find(data(:,i)>edges(end)))/ldata(i);
end

if exist('mu','var') && exist('sigma','var')
    ncol=size(data,2);
    for i=1:ncol
        dbelow(ncol+i)=normcdf(edges(  1),mu(i),sigma(i));
        dhistc(1:size(dhistc,1)-1,ncol+i)=...
            normcdf(edges(2:end  ),mu(i),sigma(i))-...
            normcdf(edges(1:end-1),mu(i),sigma(i));
        dabove(ncol+i)=norminv(edges(end),mu(i),sigma(i));
        if exist('desc','var')
            desc(ncol+i)={[desc{i} ' norm']};
        end
    end
end

%  draw the bar plot

hl1=bar(edges(1:end-1),dhistc(1:end-1,:));
ax1=gca;

%  add the annotation

title('Relative Frequency Histogram')
xlabel('Interval Edge Value')
ylabel('Relative Frequency')

if exist('desc','var')
    hleg1=legend(ax1,desc,'Location','NorthWest',...
                 'Interpreter','none');
else
    hleg1=legend(ax1);
end

%%  generate the cumulative distribution functions

% cdf=zeros(size(dhistc));
% cdf(1,:)=dhistc(1,:);
% for i=2:size(dhistc,1)
%     cdf(i,:)=cdf(i-1,:)+dhistc(i,:);
% end
cdf=cumsum(dhistc);
for i=1:size(dhistc,2)
    cdf(:,i)=dbelow(i)+cdf(:,i);
end
if exist('desc','var')
    ncol=length(desc);
    for i=1:ncol
        cdesc(i)={[desc{i} ' cdf']};
    end
end

%  draw the line plot

%  (see "Using Multiple X- and Y-Axes" and "Overlaying Other
%  Plots on Bar Graphs", or search on "YAxisLocation right")

%hold all
%hold on
%plot(edges,cdf)
%plotyy([],[],edges,cdf)

%  ticks from the bar plot will show through, so make ticks
%  from the line plot match

nytick=length(get(ax1,'YTick'));
ytmax =ceil(max(max(cdf))/0.1-0.1)*0.1;
ytinc =ytmax/(nytick-1-1);

ax2=axes('Position',get(ax1,'Position'),...
         'XLim',get(ax1,'XLim'),...
         'YLim',[0 ytinc*(nytick-1)],...
         'YTick',[0:ytinc:ytinc*(nytick-1)],...
         'XAxisLocation','bottom','YAxisLocation','right',...
         'Color','none','Layer','top');
hl2=line(edges(1:end-1),cdf(1:end-1,:),'Parent',ax2);
       
%  add the annotation

ylabel('Cumulative Percent')

% legend doesn't combine with bar chart above
if exist('cdesc','var')
    hleg2=legend(ax2,cdesc,'Location','NorthEast',...
                 'Interpreter','none');
    set(hleg2,'Color','white')
else
    hleg2=legend(ax2);
    set(hleg2,'Color','white')
end

set(gcf,'PaperPositionMode','auto')
%hold off

end
