function handles = distributionPlot(varargin)
%DISTRIBUTIONPLOT creates violin plots for convenient visualization of multiple distributions
%
% SYNOPSIS: handles = distributionPlot(data,propertyName,propertyValue,...)
%           handles = distributionPlot(ah,...)
%
% INPUT data : m-by-nData array of values, or vector of grouped data (use
%           the 'groups' property to specify the grouping variable), or
%           cell array of length nData.
%           The cell array can either contain vectors with values, or
%           m-by-2 arrays with [bins,counts] if you want to determine the
%           histograms by yourself (m can be different between cell
%           elements). Note that arrays inside cells with any
%           other shape than m-by-2 are reshaped to vector an a warning is
%           thrown (DISTRIBUTIONPLOT:AUTORESHAPE).
%
%       DISTRIBUTIONPLOT accepts the following propertyName/propertyValue
%           pairs (all are optional):
%
%       distWidth :  width of distributions; ideally between 0 and 1.
%           1 means that adjacent distributions might touch. Default: 0.9
%       variableWidth : If true, the width of the distribution changes,
%           reflecting the shape of the histogram of the data. If false,
%           the distribution is only encoded by color levels. Default: true
%       color : uniform coloring of histograms. Supply either a color
%           string ('r'), or a truecolor vector ([1 0 0]). Use a
%           cell array of length nData to specify one color per
%           distribution. Default: 'k'
%           If variableWidth is set to false, a colormap is generated that
%           goes from white to the chose color (or from black, if
%           invert==true).
%           If both 'color', and 'colormap' are specified, 'colormap' takes
%           precedence.
%       colormap : colormap used to describe the distribution (first row
%           corresponds to bins with least data, last row corresponds to
%           bins with most data (invert the grayscale colormap to have
%           black indicate the most data).
%           Supply a cell array of length nData to color distributions
%           individually. Note that using multiple colormaps means that
%           the colorbar doesn't contain much useful information.
%           Default: []
%           Colormap will index into the figure colormap, which will be
%           modified by distributionPlot. This is done to allow editing the
%           distributions in e.g. Adobe Illustrator.
%           If both 'color', and 'colormap' are specified, 'colormap' takes
%           precedence.
%       globalNorm : normalization for bin width (x-direction)
%           0 : every histogram is normalized individually so that the
%               maximum bin width is equal to distWidth. This is best
%               suited to comparing distribution shapes. Default.
%           1 : histograms are normalized such that equal bin width
%               reports equal numbers of counts per bin.
%           2 : histograms are normalized so that the relative areas
%               covered by the histograms reflect the relative total number
%               of data points.
%           3 : histograms areas are normalized so that relative densities
%               are the same across histograms. Thus, if
%               data = {rand(100,1),rand(500,1)},
%               then
%               distributionPlot(data,'globalNorm',2,'histOpt',0,'divFactor',10)
%               shows the left histogram 5x as wide as the right, while
%               distributionPlot(data,'globalNorm',3,'histOpt',0,'divFactor',10)
%               displays both histograms equally wide, since each bin
%               contains ~10% of the data.
%           Options 1 and 2 produce similar results if the bins are spaced
%           equally for the distributions. Options 0 and 3 produce similar
%           results if the data are drawn from the same distributions.
%           Note that colormaps currently always report the number of data
%           points per bin; 'globalNorm' only applies to the distribution
%           shape.
%
%       groups : grouping variable for grouped data. Grouping will be
%                   resolved by calling grp2idx, and unless xNames have
%                   been supplied, group names determine the x-labels.
%                   If the grouping variable is numeric, group labels also
%                   determine x-values, unless the parameter xValues has
%                   been specified.
%       histOpt : histogram type to plot
%                   0 : use hist command (no smoothing, fixed number of
%                       bins)
%                   1 : smoothened histogram using ksdensity with
%                       Normal kernel. Default.
%                   1.1: smoothened histogram using ksdensity where the
%                       kernel is robustly estimated via histogram.m.
%                       Normal kernel.
%                   2 : histogram command (no smoothing, automatic
%                       determination of thickness (y-direction) of bins)
%       divFactor : Parameter dependent on histOpt. If...
%                   histOpt == 0: divFactor = # of bins. Default: 25.
%                       Alternatively, pass a vector which will be
%                       interpreted as bin centers.
%                   histOpt == 1: divFactor decides by how much the default
%                       kernel-width is multiplied in order to avoid an
%                       overly smooth histogram. Default: 1/2
%                   histOpt == 2: divFactor decides by how much the
%                       automatic bin width is multiplied in order to have
%                       more (<1) or less (>1) detail. Default: 1
%       addSpread : if 1, data points are plotted with plotSpread.
%                   distWidth is ideally set to 0.95
%                   This option is not available if the data is supplied as
%                   histograms.
%                   Please download plotSpread.m separately from the File
%                   Exchange using the link in the remarks
%       addBoxes  : if 1 boxplots are overlaid on the data
%                   Default: 0
%       showMM : if 1, mean and median are shown as red crosses and
%                green squares, respectively. This is the default
%                2: only mean
%                3: only median
%                4: mean +/- standard error of the mean (no median)
%                5: mean +/- standard deviation (no median)
%                6: draw lines at the 25,50,75 percentiles (no mean)
%                0: plot neither mean nor median
%       xValues: x-coordinate where the data should be plotted.
%                If xValues are given, "distWidth" is scaled by the median
%                difference between adjacent (sorted) x-values. Note that
%                this may lead to overlapping distributions. Default:
%                1:nData
%       xNames : cell array of length nData containing x-tick names
%               (instead of the default '1,2,3')
%       xMode  : if 'auto', x-ticks are spaced automatically. If 'manual',
%                there is a tick for each distribution. If xNames is
%                provided as input, xMode is forced to 'manual'. Default:
%                'manual'.
%          NOTE: SPECIFYING XNAMES OR XVALUES OR XMODE WILL ERASE PREVIOUS
%                LABELS IF PLOTTING INTO EXISTING AXES
%       yLabel : string with label for y-axis. Default : ''
%                If empty and data is histograms, ylabel is set to 'counts'
%       invert : if 1, axes color is changed to black, and colormap is
%                   inverted.
%       histOri: Orientation of histogram. Either 'center', 'left', or
%                'right'. With 'left' or 'right', the left or right half of
%                the standard violin plot is shown. Has no effect if
%                variableWidth is false. Default: center
%       xyOri  : orientation of axes. Either 'normal' (=default), or
%                'flipped'. If 'flipped', the x-and y-axes are switched, so
%                that violin plots are horizontal. Consequently,
%                axes-specific properties, such as 'yLabel' are applied to
%                the other axis.
%       widthDiv : 1-by-2 array with [numberOfDivisions,currentDivision]
%                widthDiv allows cutting the stripe dedicated to a single
%                distribution into multible bands, which can be filled with
%                sequential calls to distributionPlot. This is one way
%                to compare two (or more) sequences of distributions. See
%                example below.
%       ah : axes handle to plot the distributions. Default: gca
%
% OUTPUT handles : 1-by-5 cell array with patch-handles for the
%                  distributions, plot handles for mean/median, the
%                  axes handle, and the plotSpread-points handle and
%                  boxplot handle
%
%
% EXAMPLES
        %--Distributions contain more information than boxplot can capture
%{
         r = rand(1000,1);
         rn = randn(1000,1)*0.38+0.5;
         rn2 = [randn(500,1)*0.1+0.27;randn(500,1)*0.1+0.73];
         rn2=min(rn2,1);rn2=max(rn2,0);
         figure
         ah(1)=subplot(3,4,1:2);
         boxplot([r,rn,rn2])
         ah(2)=subplot(3,4,3:4);
         distributionPlot([r,rn,rn2],'histOpt',2); % histOpt=2 works better for uniform distributions than the default
         set(ah,'ylim',[-1 2])

         %--- additional options  
 
         data = [randn(100,1);randn(50,1)+4;randn(25,1)+8];
         subplot(3,4,5)
 
         %--- defaults
         distributionPlot(data); 
         subplot(3,4,6)
 
         %--- show density via custom colormap only, show mean/std,
         distributionPlot(data,'colormap',copper,'showMM',5,'variableWidth',false) 
         subplot(3,4,7:8)
 
         %--- auto-binwidth depends on # of datapoints; for small n, plotting the data is useful
         % note that this option requires the additional installation
         % of plotSpread from the File Exchange (link below)
         distributionPlot({data(1:5:end),repmat(data,2,1)},'addSpread',true,'showMM',false,'histOpt',2) 
 
         %--- show quantiles
         subplot(3,4,9),distributionPlot(randn(100,1),'showMM',6)
 
         %--- horizontal orientation
         subplot(3,4,10:11),
         distributionPlot({chi2rnd(3,1000,1),chi2rnd(5,1000,1)},'xyOri','flipped','histOri','right','showMM',0),
         xlim([-3 13])
 
         %--- compare distributions side-by-side (see also example below)
         % plotting into specified axes will throw a warning that you can
         % turn off using " warning off DISTRIBUTIONPLOT:ERASINGLABELS "
         ah = subplot(3,4,12);
         subplot(3,4,12),distributionPlot(chi2rnd(3,1000,1),'histOri','right','color','r','widthDiv',[2 2],'showMM',0)
         subplot(3,4,12),distributionPlot(chi2rnd(5,1000,1),'histOri','left','color','b','widthDiv',[2 1],'showMM',0)
 
         %--Use globalNorm to generate meaningful colorbar
         data = {randn(100,1),randn(500,1)};
         figure
         distributionPlot(data,'globalNorm',true,'colormap',1-gray(64),'histOpt',0,'divFactor',[-5:0.5:5])
         colorbar

         %--Use widthDiv to compare two series of distributions
         data1 = randn(500,5);
         data2 = bsxfun(@plus,randn(500,5),0:0.1:0.4);
         figure
         distributionPlot(data1,'widthDiv',[2 1],'histOri','left','color','b','showMM',4)
         distributionPlot(gca,data2,'widthDiv',[2 2],'histOri','right','color','k','showMM',4)

         %--Christmas trees!
           x=meshgrid(1:10,1:10);
           xx = tril(x);
           xx = xx(xx>0);
           figure
           hh=distributionPlot({xx,xx,xx},'color','g','addSpread',1,'histOpt',2,'showMM',0);
           set(hh{4}{1},'color','r','marker','o')
%} 
%        END
%
% REMARKS To show distributions as clouds of points (~beeswarm plot),
%         and/or to use the option "addSpread", please download the
%         additional function plotSpread.m from the File Exchange
%         http://www.mathworks.com/matlabcentral/fileexchange/37105-plot-spread-points-beeswarm-plot
%         
%         I used to run ksdensity with the Epanechnikov kernel. However,
%         for integer data, the shape of the kernel can produce peaks
%         between the integers, which is not ideal (use histOpt=2 for
%         integer valued data).
%
%         A previous iteration of distributionPlot used the input
%         specifications below. They still work to ensure backward
%         compatibility, but are no longer supported or updated.
%           handles = distributionPlot(data,distWidth,showMM,xNames,histOpt,divFactor,invert,addSpread,globalNorm)
%           where distWidth of 1 means that the maxima
%           of  two adjacent distributions might touch. Negative numbers
%           indicate that the distributions should have constant width, i.e
%           the density is only expressed through greylevels.
%           Values between 1 and 2 are like values between 0 and 1, except
%           that densities are not expressed via graylevels. Default: 1.9
%
%
% SEE ALSO histogram, ksdensity, plotSpread, boxplot, grp2idx
%

% created with MATLAB ver.: 7.6.0.324 (R2008a) on Windows_NT
%
% created by: Jonas Dorn; jonas.dorn@gmail.com
% DATE: 08-Jul-2008
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%====================================
%% TEST INPUT
%====================================

% set defaults
def.xNames = [];
def.showMM = 1;
def.distWidth = 0.9;
def.histOpt = 1;
def.divFactor = [25,2,1];
def.invert = false;
def.colormap = [];
def.color = 'k';
def.addSpread = false;
def.addBoxes = false;
def.globalNorm = false;
def.variableWidth = true;
def.groups = [];
def.yLabel = '';
def.xValues = '';
def.xMode = 'manual';
def.histOri = 'center';
def.xyOri = 'normal';
def.widthDiv = [1 1];
isHistogram = false; %# this parameter is not set by input


if nargin == 0 || isempty(varargin{1})
    error('not enough input arguments')
end

% check for axes handle
if ~iscell(varargin{1}) && isscalar(varargin{1}) == 1 && ...
        ishandle(varargin{1}) && strcmp(get(varargin{1},'Type'),'axes')
    ah = varargin{1};
    data = varargin{2};
    varargin(1:2) = [];
    newAx = false;
    
    
else
    ah = gca;
    data = varargin{1};
    varargin(1) = [];
    newAx = true;
end

% check for current axes limits. Set NaN if the axes have no children
% yet - we need that in case we're building a complicated set of
% distributions
if ~isempty(get(ah,'children'))
    xAxLim = xlim;
    yAxLim = ylim;
else
    [xAxLim,yAxLim] = deal([NaN NaN]);
end

fh = get(ah,'Parent');

% check data. If not cell, convert
if ~iscell(data)
    [nPoints,nData] = size(data);
    data = mat2cell(data,nPoints,ones(nData,1));
else
    % get nData
    data = data(:);
    nData = length(data);
    % make sure all are vectors
    badCol = ~cellfun(@isvector,data) & ~cellfun(@isempty,data);
    if any(badCol)
        nCols = cellfun(@(x)(size(x,2)),data(badCol));
        if all(nCols==2)
            % bins,counts
            isHistogram = true;
        else
            warning('DISTRIBUTIONPLOT:AUTORESHAPE',...
                'Elements %s of the cell array are not vectors. They will be reshaped automatically',...
                num2str(find(badCol)'));
            data(badCol) = cellfun(@(x)(x(:)),data(badCol),'UniformOutput',false);
        end
    end
end

parserObj = inputParser;
parserObj.FunctionName = 'distributionPlot';
stdWidth = 1; % scaling parameter for variableWidth with uneven x-values
% check whether we're dealing with pN/pV or straight arguments
if ~isempty(varargin) && ~ischar(varargin{1}) && ~isstruct(varargin{1})
    % use old format
    % distWidth,showMM,xNames,histOpt,divFactor,invert,addSpread,globalNorm
    def.distWidth = 1.9;
    parserObj.addOptional('distWidth',def.distWidth);
    parserObj.addOptional('showMM',def.showMM);
    parserObj.addOptional('xNames',def.xNames);
    parserObj.addOptional('histOpt',def.histOpt);
    parserObj.addOptional('divFactor',def.divFactor);
    parserObj.addOptional('invert',def.invert);
    parserObj.addOptional('addSpread',def.addSpread);
    parserObj.addOptional('globalNorm',def.globalNorm);
    parserObj.addOptional('groups',def.groups);
    parserObj.addOptional('yLabel',def.yLabel);
    parserObj.addOptional('color',def.color);
    
    
    parserObj.parse(varargin{:});
    opt = parserObj.Results;
    % fill in defaults that are not supported in the old version of the
    % code
    opt.colormap = [];
    opt.variableWidth = true;
    opt.histOri = 'center';
    opt.xValues = [];
    opt.xMode = 'auto';
    opt.xyOri = 'normal';
    opt.widthDiv = [1 1];
    
    % overwrite empties with defaults - inputParser considers empty to be a
    % valid input.
    fnList = fieldnames(opt);
    for fn = fnList'
        if isempty(opt.(fn{1}))
            opt.(fn{1}) = def.(fn{1});
        end
    end
    
    
    % fix a few parameters
    if opt.distWidth > 1
        opt.distWidth = opt.distWidth - 1;
    else
        opt.colormap = 1-gray(128);
    end
    if opt.distWidth < 0
        opt.variableWidth = false;
        opt.distWidth = abs(opt.distWidth);
    end
    
    if ~isempty(opt.xNames)
        opt.xMode = 'manual';
    end
    
    
else
    defNames = fieldnames(def);
    for dn = defNames(:)'
        try
        parserObj.addParamValue(dn{1},def.(dn{1})); %#ok<NVREPL>
        catch
            parserObj.addParameter(dn{1},def.(dn{1}));
        end
    end
    
    
    parserObj.parse(varargin{:});
    opt = parserObj.Results;
    
    % if groups: deal with data
    if ~isempty(opt.groups)
        [idx,labels,vals] = grp2idx(opt.groups);
        % convert data to cell array
        data = accumarray(idx,data{1},[],@(x){x});
        nData = length(data);
        % if not otherwise provided, use group labels for xnames
        if isempty(opt.xNames)
            opt.xNames = labels;
            if ~iscell(opt.xNames)
                opt.xNames = num2cell(opt.xNames);
            end
        end
        if isnumeric(vals) && isempty(opt.xValues)
            opt.xValues = vals;
        end
        
    end
    
    if ~ischar(opt.xyOri) || ~any(ismember(opt.xyOri,{'normal','flipped'}))
        error('option xyOri must be either ''normal'' or ''flipped'' (is ''%s'')',opt.xyOri);
    end
    
    
    
    
end
% common checks

% default x-values: 1:n
if isempty(opt.xValues)
    opt.xValues = 1:nData;
elseif length(opt.xValues) ~= nData
    error('please supply as many x-data values as there are data entries')
elseif length(opt.xValues) > 1 % only check for scale if more than 1 value
    % scale width
    stdWidth = median(diff(sort(opt.xValues)));
    opt.distWidth = opt.distWidth * stdWidth;
end


if ~isscalar(opt.divFactor) && length(opt.divFactor) == 3 && all(opt.divFactor==def.divFactor)
    opt.divFactor = opt.divFactor(floor(opt.histOpt)+1);
end
if isHistogram
    opt.histOpt = 99;
    if isempty(opt.yLabel)
        opt.yLabel = 'counts';
    end
end



% check colors/colormaps: do we need to expand colormap?
if ~iscell(opt.colormap)
    opt.colormap = {opt.colormap};
end
if ~iscell(opt.color)
    opt.color = {opt.color};
end
for iColor = 1:length(opt.color)
    if ischar(opt.color{iColor})
        opt.color{iColor} = colorCode2rgb(opt.color{iColor});
    end
end

% expand - if only single colormap specified, we expand only once
if ~opt.variableWidth
    missingColormaps = find(cellfun(@isempty,opt.colormap));
    for iMissing = missingColormaps(:)'
        
        endColor = opt.color{max(iMissing,length(opt.color))};
        % normally, we go from white to color
        cmap = zeros(128,3);
        for rgb = 1:3
            cmap(:,rgb) = linspace(1,endColor(rgb),128);
        end
        opt.colormap{iMissing} = cmap;
        
    end
end

% if we have colormaps, we need to create a master which we add to the
% figure. Invert if necessary, and expand the cell array to nData
colormapLength = cellfun(@(x)size(x,1),opt.colormap);
if any(colormapLength>0)
    
    colormap = cat(1,opt.colormap{:});
    if opt.invert
        colormap = 1-colormap;
    end
    set(fh,'Colormap',colormap)
    if length(opt.colormap) == 1
        opt.colormap = repmat(opt.colormap,nData,1);
        colormapLength = repmat(colormapLength,nData,1);
        colormapOffset = zeros(nData,1);
        singleMap = true;
    else
        colormapOffset = [0;cumsum(colormapLength(1:end-1))];
        singleMap = false;
    end
    
else
    
    colormapLength = zeros(nData,1);
    if length(opt.color) == 1
        opt.color = repmat(opt.color,nData,1);
    end
    if opt.invert
        opt.color = cellfun(@(x)1-x,opt.color,'uniformOutput',false);
    end
end


% set hold on
holdState = get(ah,'NextPlot');
set(ah,'NextPlot','add');

% if new axes: invert
if newAx && opt.invert
    set(ah,'Color','k')
end

%===================================



%===================================
%% PLOT DISTRIBUTIONS
%===================================

% assign output
hh = NaN(nData,1);
[m,md,sem,sd] = deal(nan(nData,1));
if opt.showMM == 6
    md = nan(nData,3,3); % md/q1/q3, third dim is y/xmin/xmax
end

% make sure xValues are not something weird, like an enum
opt.xValues = double(opt.xValues);


% get base x-array
% widthDiv is a 1-by-2 array with
% #ofDivs, whichDiv
% The full width (distWidth) is split into
% #ofDivs; whichDiv says which "stripe" is active
xWidth = opt.distWidth/opt.widthDiv(1);
xMin = -opt.distWidth/2;
xLow = xMin + xWidth * (opt.widthDiv(2)-1);
xBase = [-xWidth;xWidth;xWidth;-xWidth]/2;
xOffset = xLow + xWidth/2;

% b/c of global norm: loop twice
plotData = cell(nData,2);

% loop through data. Prepare patch input, then draw patch into gca
for iData = 1:nData
    currentData = data{iData};
    % only plot if there is some finite data
    if ~isempty(currentData(:)) && any(isfinite(currentData(:)))
        
        switch floor(opt.histOpt)
            case 0
                % use hist
                [xHist,yHist] = hist(currentData,opt.divFactor);
                
            case 1
                % use ksdensity
                
                if opt.histOpt == 1.1
                    % use histogram to estimate kernel
                    [dummy,x] = myHistogram(currentData); %#ok<ASGLU>
                    if length(x) == 1
                        % only one value. Make fixed distribution
                        dx = 0.1;
                        yHist = x;
                        xHist = sum(isfinite(currentData));
                    else
                        dx = x(2) - x(1);
                    
                    % make sure we sample frequently enough
                    x = min(x)-dx:dx/3:max(x)+dx;
                    [xHist,yHist] = ksdensity(currentData,x,'kernel','normal','width',dx/(1.5*opt.divFactor));
                    end
                else
                    
                    % x,y are switched relative to normal histogram
                    [xHist,yHist,u] = ksdensity(currentData,'kernel','normal');
                    % take smaller kernel to avoid over-smoothing
                    if opt.divFactor ~= 1
                        [xHist,yHist] = ksdensity(currentData,'kernel','normal','width',u/opt.divFactor);
                    end
                end
                
                % modify histogram such that the sum of bins (not the
                % integral under the curve!) equals the total number of
                % observations, in order to be comparable to hist
                xHist = xHist/sum(xHist)*sum(isfinite(currentData));
                
            case 2
                % use histogram - bar heights are counts as in hist
                [xHist,yHist] = myHistogram(currentData,opt.divFactor,0);
            case 99
                % bins,counts already supplied
                xHist = currentData(:,2)';
                yHist = currentData(:,1)';
        end
        plotData{iData,1} = xHist;
        plotData{iData,2} = yHist;
    end
end

goodData = find(~cellfun(@isempty,plotData(:,1)));
% get norm
switch opt.globalNorm
    case 3
        % #3 normalizes relative densities
        xNorm(goodData) = cellfun(@(x)min(diff(x)),plotData(goodData,2));
        xNorm(goodData) = xNorm(goodData) .* cellfun(@sum,plotData(goodData,1))';
        maxNorm(goodData) = cellfun(@max,plotData(goodData,1));
        xNorm(goodData) = xNorm(goodData)*max(maxNorm(goodData)./xNorm(goodData));
        
    case 2
        % #2 should normalize so that the integral of the
        % different histograms (i.e. area covered) scale with the
        % respective sum of counts across all bins. Requires evenly spaced
        % histograms at the moment
        xNorm(goodData) = cellfun(@(x)min(diff(x)),plotData(goodData,2));
        maxNorm(goodData) = cellfun(@max,plotData(goodData,1));
        xNorm(goodData) = xNorm(goodData)*max(maxNorm(goodData)./xNorm(goodData));
    case 1
        xNorm(goodData) = max(cat(2,plotData{:,1}));
    case 0
        xNorm(goodData) = cellfun(@max,plotData(goodData,1));
end


for iData = goodData'
    
    % find current data again
    currentData = data{iData};
    
    xHist = plotData{iData,1};
    yHist = plotData{iData,2};
    
    % find y-step
    dy = min(diff(yHist));
    if isempty(dy)
        dy = 0.1;
    end
    
    % create x,y arrays
    nPoints = length(xHist);
    xArray = repmat(xBase,1,nPoints);
    yArray = repmat([-0.5;-0.5;0.5;0.5],1,nPoints);
    
    
    % x is iData +/- almost 0.5, multiplied with the height of the
    % histogram
    if opt.variableWidth
        
        
        tmp = xArray.*repmat(xHist,4,1)./xNorm(iData);
        
        switch opt.histOri
            case 'center'
                % we can simply use xArray
                xArray = tmp;
            case 'right'
                % shift everything to the left
                delta = tmp(1,:) - xArray(1,:);
                xArray = bsxfun(@minus,tmp,delta);
            case 'left'
                % shift everything to the right
                delta = tmp(1,:) - xArray(1,:);
                xArray = bsxfun(@plus,tmp,delta);
        end
        
        xArray = xArray + opt.xValues(iData);
        
    else
        xArray = xArray + iData;
    end
    
    % add offset (in case we have multiple widthDiv)
    xArray = xArray + xOffset;
    
    
    % yData is simply the bin locations
    yArray = repmat(yHist,4,1) + dy*yArray;
    
    % add patch
    vertices = [xArray(:),yArray(:)];
    faces = reshape(1:numel(yArray),4,[])';
    
    if colormapLength(iData) == 0
        colorOpt = {'FaceColor',opt.color{iData}};
    else
        % calculate index into colormap
        if singleMap
            % use scaled mapping so that colorbar is meaningful
            if opt.globalNorm > 0
                colorOpt = {'FaceVertexCData',xHist','CDataMapping','scaled','FaceColor','flat'};
            else
                colorOpt = {'FaceVertexCData',xHist'/xNorm(iData),'CDataMapping','scaled','FaceColor','flat'};
            end
            
        else
            idx = round((xHist/xNorm(iData))*(colormapLength(iData)-1))+1;
            colorOpt = {'FaceVertexCData',idx'+colormapOffset(iData),'CDataMapping','direct','FaceColor','flat'};
        end
    end
    
    
    switch opt.xyOri
        case 'normal'
            hh(iData)= patch('Vertices',vertices,'Faces',faces,'Parent',ah,colorOpt{:},'EdgeColor','none');
        case 'flipped'
            hh(iData)= patch('Vertices',vertices(:,[2,1]),'Faces',faces,'Parent',ah,colorOpt{:},'EdgeColor','none');
    end
    
    if opt.showMM > 0
        if isHistogram
            [m(iData),sem(iData)] = weightedStats(currentData(:,1),currentData(:,2),'w');
            sd(iData) = sem(iData) * sqrt(sum(currentData(:,2)));
            % weighted median: where we're at middle weight
            % may need some tweaking
            goodCurrentData = sortrows(currentData(all(isfinite(currentData),2),:),1);
            weightList = cumsum(goodCurrentData(:,2));
            weightList = weightList / weightList(end);
            md(iData) = goodCurrentData(find(weightList>0.5,1,'first'),1);
        else
            m(iData) = nanmean(currentData);
            md(iData) = nanmedian(currentData);
            sd(iData) = nanstd(currentData);
            sem(iData) = sd(iData)/sqrt(sum(isfinite(currentData)));
        end
        
        if opt.showMM == 6
            % read quantiles - "y"-value, plus x-start-stop
            % re-use md array which allows using a loop below instead of
            % lots of copy-paste
            % md array is md/q1/q3, with third dimension y/xmin/xmax
            
            md(iData,2,1) = prctile(currentData,25);
            md(iData,3,1) = prctile(currentData,75);
            
            for qq = 1:3
                % find corresponding y-bin
                yLoc =  repmat(...
                    any(yArray>md(iData,qq,1),1) & any(yArray<=md(iData,qq,1),1),...
                    [4 1]);
                % look up corresponding x-values. Note that there is a bit
                % of a risk that the line will be exactly between two very
                % different bins - but if we make the line longer, it will
                % be ugly almost all the time
                md(iData,qq,2) = min( xArray( yLoc ) );
                md(iData,qq,3) = max( xArray( yLoc ) );
            end
            
        end
    end
end % loop

sh = [];
if opt.addSpread
    if isHistogram
        disp('Option addSpread is unavailable if data is supplied as histograms. Call plotSpread separately')
    else
        % add spread
        try
            sh = plotSpread(ah,data,'xValues',opt.xValues,'xyOri',opt.xyOri);
            set(sh{1}(ishandle(sh{1})),'color',[0,128,255]/255);
        catch me
            if strcmp(me.identifier,'MATLAB:UndefinedFunction')
                error('plotSpread not found. Please download it from the Matlab File Exchange')
            else
                rethrow(me)
            end
        end
    end
end


bh = [];
if opt.addBoxes
    if isHistogram
        disp('Option addBox is unavailable if data is supplied as histograms. Call boxPlot separately')
    else
        % add spread
        try
            switch opt.xyOri
                case 'normal'
                    ori = 'vertical';
                case 'flipped'
                    ori = 'horizontal';
            end
            % boxplot wants a matrix. fix up data
            matHeight = max(cellfun(@length,data));
            boxMat = NaN(matHeight,nData);
            for iData = 1:nData
                boxMat(1:length(data{iData}),iData) = data{iData}(:);
            end
            bh = boxplot(ah,boxMat,'labels',opt.xValues,'orientation',ori);
        catch me
            if strcmp(me.identifier,'MATLAB:UndefinedFunction')
                error('plotSpread not found. Please download it from the Matlab File Exchange')
            else
                rethrow(me)
            end
        end
    end
end



mh = [];mdh=[];
if opt.showMM
    % plot mean, median. Mean is filled red circle, median is green square
    % I don't know of a very clever way to flip xy and keep everything
    % readable, thus it'll be copy-paste
    switch opt.xyOri
        case 'normal'
            if any(opt.showMM==[1,2])
                mh = plot(ah,opt.xValues+xOffset,m,'+r','Color','r','MarkerSize',12);
            end
            if any(opt.showMM==[1,3])
                mdh = plot(ah,opt.xValues+xOffset,md,'sg','MarkerSize',12);
            end
            if opt.showMM == 4
                mh = plot(ah,opt.xValues+xOffset,m,'+r','Color','r','MarkerSize',12);
                mdh = myErrorbar(ah,opt.xValues+xOffset,m,sem);
            end
            if opt.showMM == 5
                mh = plot(ah,opt.xValues+xOffset,m,'+r','Color','r','MarkerSize',12);
                mdh = myErrorbar(ah,opt.xValues+xOffset,m,sd);
            end
            if opt.showMM == 6
                mdh(1,:) = plot(ah,squeeze(md(:,1,2:3))',repmat(md(:,1,1)',2,1),'color','r','lineWidth',2);%,'lineStyle','--');
                mdh(2,:) = plot(ah,squeeze(md(:,2,2:3))',repmat(md(:,2,1)',2,1),'color','r','lineWidth',1);%,'lineStyle','--');
                mdh(3,:) = plot(ah,squeeze(md(:,3,2:3))',repmat(md(:,3,1)',2,1),'color','r','lineWidth',1);%,'lineStyle','--');
            end
        case 'flipped'
            if any(opt.showMM==[1,2])
                mh = plot(ah,m,opt.xValues+xOffset,'+r','Color','r','MarkerSize',12);
            end
            if any(opt.showMM==[1,3])
                mdh = plot(ah,md,opt.xValues+xOffset,'sg','MarkerSize',12);
            end
            if opt.showMM == 4
                mh = plot(ah,m,opt.xValues+xOffset,'+r','Color','r','MarkerSize',12);
                mdh = myErrorbar(ah,m,opt.xValues+xOffset,[sem,NaN(size(sem))]);
            end
            if opt.showMM == 5
                mh = plot(ah,m,opt.xValues+xOffset,'+r','Color','r','MarkerSize',12);
                mdh = myErrorbar(ah,m,opt.xValues+xOffset,[sd,NaN(size(sd))]);
            end
            if opt.showMM == 6
                mdh(1,:) = plot(ah,repmat(md(:,1,1)',2,1),squeeze(md(:,1,2:3))','color','r','lineWidth',2);%,'lineStyle','--');
                mdh(2,:) = plot(ah,repmat(md(:,2,1)',2,1),squeeze(md(:,2,2:3))','color','r','lineWidth',1);%,'lineStyle','--');
                mdh(3,:) = plot(ah,repmat(md(:,3,1)',2,1),squeeze(md(:,3,2:3))','color','r','lineWidth',1);%,'lineStyle','--');
            end
    end
end

% find extents of x-axis (or y-axis, if flipped)
minX = min(opt.xValues)-stdWidth;
maxX = max(opt.xValues)+stdWidth;

if ~isnan(xAxLim(1))
    % we have previous limits
    switch opt.xyOri
        case 'normal'
            minX = min(minX,xAxLim(1));
            maxX = max(maxX,xAxLim(2));
        case 'flipped'
            minX = min(minX,yAxLim(1));
            maxX = max(maxX,yAxLim(2));
    end
end


% if ~empty, use xNames
switch opt.xyOri
    case 'normal'
        switch opt.xMode
            case 'manual'
                if newAx == false
                    warning('DISTRIBUTIONPLOT:ERASINGLABELS','Plotting into an existing axes and specifying labels will erase previous labels')
                end
                % x-values are sorted when plotting
                [sortedX,sortIdx] = sort(opt.xValues);
                set(ah,'XTick',sortedX);
                if ~isempty(opt.xNames)
                    set(ah,'XTickLabel',opt.xNames(sortIdx))
                end
            case 'auto'
                % reset to auto
                set(ah,'XTickMode','auto')
        end
        if ~isempty(opt.yLabel)
            ylabel(ah,opt.yLabel);
        end
        % have plot start/end properly
        xlim([minX,maxX])
    case 'flipped'
        switch opt.xMode
            case 'manual'
                if newAx == false
                    warning('DISTRIBUTIONPLOT:ERASINGLABELS','Plotting into an existing axes and specifying labels will erase previous labels')
                end
               % x-values are sorted when plotting
                [sortedX,sortIdx] = sort(opt.xValues);
                set(ah,'YTick',sortedX);
                if ~isempty(opt.xNames)
                    set(ah,'YTickLabel',opt.xNames(sortIdx))
                end
            case 'auto'
                % no need to do anything
        end
        if ~isempty(opt.yLabel)
            xlabel(ah,opt.yLabel);
        end
        % have plot start/end properly
        ylim([minX,maxX])
end


%==========================


%==========================
%% CLEANUP & ASSIGN OUTPUT
%==========================

if nargout > 0
    handles{1} = hh;
    handles{2} = [mh;mdh];
    handles{3} = ah;
    handles{4} = sh;
    handles{5} = bh;
end

set(ah,'NextPlot',holdState);