function varargout = beeswarm_plot(varargin)

% beeswarm_plot: Scatter plot of data with random jitter added to the data
% for better visualisation. A notched boxplot is superimposed on the
% scatter plot. A t-test is done between adjoining groups and significant
% p-value is indicated with asterisk. 
% 
% beeswarm_plot({data1}) makes a beeswarm plot of the data in data1.
% 
% beeswarm_plot({data1 data2}) makes a beeswarm plot of data1 & data2.
% 
% beeswarm_plot({data1 data2; data3 data4}) makes error plot of data1 & 
% data2 together and data3 & data4 together
% 
% beeswarm_plot(...,'color',[a b c]) makes a beeswarm_plot using the color
% mentioned in the color vector
% 
% beeswarm_plot(...,'stats','x') plots the descriptive stats of the data.
% If the string is 'mean' then mean+-sem is drawn. If the string is
% 'boxplot' then a notched box-plot is drawn. the default is box-plot.
% 
% beeswarm_plot(...,'scatter','x') plots the raw data points.
% If the string is 'random' then the raw data points are randomly
% scattered. If the string is 'lined' then the raw data points are plotted
% as a histogram.
% 
% beeswarm_plot(....'distribution','x') scatters the raw data points if
% 'scatter' is set as 'lined'. 
% If the string is 'hist' then the raw data points are randomly scattered
% in each bin of the histogram. If the string is 'leaf' then the data
% points are plotted like a leaf with points having the smallest or largest
% deviation from the median of the mean having the least shift.
% 
% beeswarm_plot(...,'tail',tail) does a two-tailed test if 'both' or does
% one-sided test if 'right' or 'left'
% 
% Input variables
% ---------------
% The data to be plotted as an cell array
% 
% 
% Optional input arguments
% ------------------------
% color: RGB values of the colors
% 
% plottype: whether 'mean' or 'boxplot'
% 
% scatter: whether 'random' or 'lined'
% 
% distribution: whether 'hist' or 'leaf'
% 
% tail: whether 'left', 'right', or 'both'. Default = 'both'
% 
% stats: whether 'param' or 'nonparam'. Default = 'default'
% 
% ********************************************************************
% Conjured by Sumitash
% 
% 5 days after a splendid Valentine's Day, 2014
% Added matrix plotting 2 weeks after India's Independence Day, 2018
% Fixed CI median on the day the Nobel Prize was estabilished, 2018
% ********************************************************************

% --- Input ---
val_inp = varargin{1,1};

% --- Stats test type ---
istest = cellfun(@(x) ischar(x) && strcmp(x, 'stats'), varargin);
if any(istest)
    varargin = varargin(~istest);
    usetest = varargin{istest};
    if strcmp(usetest,'param')==0 & strcmp(usetest,'nonparam')==0
        warning('U non-normal moron! Not recognisable test type! Using default test! B-)');
        usetest = 'default';
    end
else
    usetest = 'default';
end 

% --- Tails of the stat test ---
istail = cellfun(@(x) ischar(x) && strcmp(x, 'tail'), varargin);
if any(istail)
    varargin = varargin(~istail);
    usetail = varargin{istail};
else
    usetail = 'both';
end

% --- Stats type ---
isplot = cellfun(@(x) ischar(x) && strcmp(x, 'plottype'), varargin);
if any(isplot)
    varargin = varargin(~isplot);
    useplot = cell2mat((varargin(isplot)));
    if strcmp(useplot,'boxplot')==0 & strcmp(useplot,'mean')==0
        warning('U doofus! Not recognisable stats method! Using boxplot! B-)');
        useplot = 'boxplot';
    end
else
    useplot = 'boxplot';
end

% --- Scatter type ---
isscatter = cellfun(@(x) ischar(x) && strcmp(x, 'scatter'), varargin);
if any(isscatter)
    varargin = varargin(~isscatter);
    usescatter = varargin{isscatter};
    if strcmp(usescatter,'random')==0 & strcmp(usescatter,'lined')==0
        warning('U scatter-brain! Not recognisable scatter type! Using lined points! B-)');
        usescatter = 'lined';
    end
else
    usescatter = 'lined';
end 

% --- Distribution type ---
isdistb = cellfun(@(x) ischar(x) && strcmp(x, 'distribution'), varargin);
if any(isdistb)
    varargin = varargin(~isdistb);
    usedistb = varargin{isdistb};
    if strcmp(usedistb,'hist')==0 & strcmp(usedistb,'leaf')==0
        warning('U block-head! Not recognisable distribution type! Using hist! B-)');
        usedistb = 'hist';
    end
else
    usedistb = 'hist';
end 
% Extract values
if size(val_inp,1)==1
    for i = 1:size(val_inp,2)
        size_val(i) = length(cell2mat(val_inp(1,i)));
    end
    val = NaN(size(val_inp,2),max(size_val));
    for i = 1:size(val_inp,2)
        val(i,1:length(val_inp{1,i})) = val_inp{1,i};
    end
    plot_matrix = 0;
else
    count       = 1;
    size_val    = max(cellfun('length',val_inp),[],2);
    val         = NaN(size(val_inp,1)*size(val_inp,2),max(size_val));
    for row = 1:size(val_inp,1)
        for col = 1:size(val_inp,2)
            eval(['val(' num2str(count) ',1:length(val_inp{row,col})) = val_inp{row,col};'])
            count = count + 1;
        end
    end
    plot_matrix = 1;
end
max_val = max(size(val));

% --- Color ---
iscolor = cellfun(@(x) ischar(x) && strcmp(x, 'color'), varargin);
if any(iscolor)
    varargin = varargin(~iscolor);
    usecolor = cell2mat((varargin(iscolor)));
else
    usecolor = lines(size(val,1));
end

% --- Marker size (default - 6) ---
issize = cellfun(@(x) ischar(x) && strcmp(x, 'markersize'), varargin);
if any(issize)
    varargin = varargin(~issize);
    usesize  = cell2mat(varargin(issize));
    if ~isnumeric(usesize)
        warning('U missized fool! Forgot your size! Using size 6! B-)');
        usesize = 40;
    end
else
    usesize  = 40;
end

% --- Calculate histogram ---
no_of_bins = 10;
for i = 1:size(val,1)
    edges = linspace(min(min(val)),max(max(val)),no_of_bins);
    [n(i,:),binindx(i,:)] = histc(val(i,:),edges);
end
while any(any(n>50)) & no_of_bins<40                                       % increase the number of bins if n>40
    no_of_bins = no_of_bins + 1;
    edges = linspace(min(min(val)),max(max(val)),no_of_bins);
    n = []; binindx = [];
    for i = 1:size(val,1)
        [n(i,:),binindx(i,:)] = histc(val(i,:),edges);
    end
end

% --- Find scattering position ---
if strcmp(usescatter,'lined')
    scat_dif = 0.015;
    if (plot_matrix==0) 
        if size(val_inp,2)>2 
            xscatter = 0.015;
        else
            xscatter = 0.02;
        end
    else
        if size(val_inp,2)>2 
            xscatter = 0.004;
        else
            xscatter = 0.006;
        end
    end
    scat_val = NaN(size(val));
    for i = 1:size(val,1) 
        for j = 1:size(n,2)
            if n(i,j)>0
                indx = find(binindx(i,:)==j);
                if strcmp(usedistb,'leaf')==1
                    dummy = sortrows([((val(i,indx)-nanmedian(val(i,indx))))' indx']);                % sort rt depending on deviation from median of the bin
                    indx = dummy(:,2)';
                    dum = -scat_dif*floor(length(indx)/2):scat_dif:scat_dif*floor(length(indx)/2);
                    dum = sort(abs(dum));
                    dum(2:2:end) = dum(2:2:end).*-1;
                    if mod(length(indx),2)==0                                                       % remove 0 to fight the demon of dimension mismatch
                        dum = dum(dum~=0);
                    end
                    scat_val(i,indx) = dum;
                else
                    scat_val(i,indx) = (xscatter:xscatter:n(i,j)*xscatter) - ((n(i,j)+1)*xscatter)/2;
                end
            end
        end
    end
end

% --- Find scattering factor ---
if strcmp(usescatter,'random')
    if max(size_val)<50
        scat_factor = 30;
    elseif max(size_val)>=50 & max(size_val)<100
        scat_factor = 20;
    else
        scat_factor = 10;
    end
end

% --- Plot scatter plot ---
if (plot_matrix==0) 
    xpos = 1:size(val_inp,2);
    if size(val_inp,2)>2 
        boxwidth = 0.15;
        medwidth = 0.1;
    else
        boxwidth = 0.2;
        medwidth = 0.15;
    end
else
    xpos = ones(size(val_inp));
    xpos = bsxfun(@times,xpos,(1:size(val_inp,1))');
    if size(val_inp,2)==2
        shift_xpos = [-0.15 0.15];
    elseif size(val_inp,2)==3
        shift_xpos = [-0.225 0 0.225];
    elseif size(val_inp,2)==4
        shift_xpos = [-0.3 -0.1 0.1 0.3];
    elseif size(val_inp,2)==5
        shift_xpos = [-0.4 -0.2 0 0.2 0.4];
    end
    xpos = bsxfun(@plus,xpos,shift_xpos)';
    xpos = xpos(:);

    if size(val_inp,2)>2 
        boxwidth = 0.08;
        medwidth = 0.05;
    else
        boxwidth = 0.1;
        medwidth = 0.06;
    end
end
for i = 1:size(val,1)
    if strcmp(usescatter,'random')
        scat_val = (repmat(xpos(i),1,length(val(i,:)))+randn(1,length(val(i,:))))./scat_factor;       
        h(i) = scatter(scat_val,val(i,:),usesize,usecolor(i,1:3),'filled','MarkerEdgeColor','k');    % make scatter plot
    else
        scat_val(i,:) = scat_val(i,:) + xpos(i);
        h(i) = scatter(scat_val(i,:),val(i,:),usesize,usecolor(i,1:3),'filled','MarkerEdgeColor','k');
    end
    hold on
    if strcmp(useplot,'mean')
        plot([xpos(i)-boxwidth xpos(i)+boxwidth],repmat(nanmean(val(i,:)),1,2),'k-','LineWidth',2);  % mean
        plot([xpos(i)-boxwidth xpos(i)+boxwidth],repmat(nanmean(val(i,:))-nanstd(val(i,:),1)./sqrt(sum(~isnan(val(i,:)))),1,2),'k-','LineWidth',1.25); % sem
        plot([xpos(i)-boxwidth xpos(i)+boxwidth],repmat(nanmean(val(i,:))+nanstd(val(i,:),1)./sqrt(sum(~isnan(val(i,:)))),1,2),'k-','LineWidth',1.25); % sem
    elseif strcmp(useplot,'boxplot')
        % --- stats for box-plot ---
        prctile_rt = prctile(val(i,:),[25 50 75]);
        ci_median = [prctile_rt(2)-1.7*(1.25*(prctile_rt(3)-prctile_rt(1))/(1.35*sqrt(sum(~isnan(val(i,:)))))) ...
                     prctile_rt(2)+1.7*(1.25*(prctile_rt(3)-prctile_rt(1))/(1.35*sqrt(sum(~isnan(val(i,:))))))];
        lowlim = max((prctile_rt(1)-1.5*iqr(val(i,:))),nanmin(val(i,:)));             % lower limit for outliers
        upplim = min((prctile_rt(3)+1.5*iqr(val(i,:))),nanmax(val(i,:)));             % upper limit for outliers
        
        plot([xpos(i)-medwidth xpos(i)+medwidth],repmat(prctile_rt(2),1,2),'k-','LineWidth',3);         % median
        
        % --- CI of median ---
        plot([xpos(i)-boxwidth xpos(i)-medwidth],[ci_median(1) prctile_rt(2)],'k-','LineWidth',1.5);     
        plot([xpos(i)-boxwidth xpos(i)-medwidth],[ci_median(2) prctile_rt(2)],'k-','LineWidth',1.5);
        plot([xpos(i)+boxwidth xpos(i)+medwidth],[ci_median(1) prctile_rt(2)],'k-','LineWidth',1.5);
        plot([xpos(i)+boxwidth xpos(i)+medwidth],[ci_median(2) prctile_rt(2)],'k-','LineWidth',1.5);
        % --- box-plot lines ---
        plot([xpos(i)-boxwidth xpos(i)+boxwidth],repmat(prctile_rt(1),1,2),'k-','LineWidth',1.5);
        plot([xpos(i)-boxwidth xpos(i)+boxwidth],repmat(prctile_rt(3),1,2),'k-','LineWidth',1.5);
        plot([xpos(i)-boxwidth xpos(i)-boxwidth],[ci_median(2) prctile_rt(3)],'k-','LineWidth',1.5);
        plot([xpos(i)+boxwidth xpos(i)+boxwidth],[ci_median(2) prctile_rt(3)],'k-','LineWidth',1.5);
        plot([xpos(i)-boxwidth xpos(i)-boxwidth],[ci_median(1) prctile_rt(1)],'k-','LineWidth',1.5);
        plot([xpos(i)+boxwidth xpos(i)+boxwidth],[ci_median(1) prctile_rt(1)],'k-','LineWidth',1.5);
        % --- errorbars of box-plot ---
        errorbar(xpos(i),prctile_rt(3),0,upplim-prctile_rt(3),'k-','LineWidth',1.5)
        errorbar(xpos(i),prctile_rt(1),lowlim-prctile_rt(1),0,'k-','LineWidth',1.5)
    end
end

% --- Put p-value significance marks ---
axes_pos = get(gca,'Position');
if (plot_matrix==0)
    if size(val_inp,2)>2 
        start_arrow = axes_pos(1) + (axes_pos(3)/(size(val,1)-1))*(0.3 + (size(val,1))*0.05);     % Starting point of double arrow
        add_arrow = axes_pos(1) + (axes_pos(3)/(size(val,1)-1))*(0.7 + (size(val,1))*0.05);       % Ending point of double arrow
    else
        start_arrow = axes_pos(1) + (axes_pos(3)/(size(val,1)-1))*0.32;     % Starting point of double arrow
        add_arrow = axes_pos(1) + (axes_pos(3)/(size(val,1)-1))*0.68;       % Ending point of double arrow
    end
else
%     if size(val_inp,2)>2
%         start_arrow = axes_pos(1) + (axes_pos(3)/(size(val,1)-1))*(0.07 + ((size(val,1))-1)*0.04);     % Starting point of double arrow
%         add_arrow = axes_pos(1) + (axes_pos(3)/(size(val,1)-1))*(0.20 + ((size(val,1))-1)*0.04);       % Ending point of double arrow
%     else
        start_arrow = axes_pos(1) + (axes_pos(3)/(size(val,1)-1))*0.6;     % Starting point of double arrow
        add_arrow = axes_pos(1) + (axes_pos(3)/(size(val,1)-1))*0.9;       % Ending point of double arrow
%     end
end

if (plot_matrix==0)
    for i = 1:size(val,1)-1
        if strcmp(usetest,'default')
            if (sum(~isnan(val(i,:)))>4 & sum(~isnan(val(i+1,:)))>4)
                if ((lillietest(val(i,:))==0 & lillietest(val(i+1,:))==0 | ...
                        sum(~isnan(val(i,:)))>30 & sum(~isnan(val(i+1,:)))>30))
                    [~,p(i)] = ttest2(val(i,:),val(i+1,:),0.05,usetail);
                else
                    p(i) = ranksum(val(i,:),val(i+1,:));
                end
            else
                try
                    p(i) = ranksum(val(i,:),val(i+1,:));
                catch
                    p(i) = NaN;
                end
            end
        elseif strcmp(usetest,'param')
            [~,p(i)] = ttest2(val(i,:),val(i+1,:),0.05,usetail);
        elseif strcmp(usetest,'nonparam')
            p(i) = ranksum(val(i,:),val(i+1,:));
        end

        if p(i) <= 0.05 & p(i) >= 0.01
            annotation('textbox',[mean([start_arrow add_arrow])-0.02 axes_pos(2)/3+axes_pos(4) 0.1 0.05],'String','*','FontSize',20,'FontWeight','demi','EdgeColor','none');
%             annotation('textbox',[mean([start_arrow add_arrow])-0.09 axes_pos(2)/3+axes_pos(4) 0.2 0.05],'String',strcat('p =',sprintf(' %0.3f',p(i))),'FontSize',10,'FontWeight','demi','EdgeColor','none');
        elseif p(i) < 0.01 & p(i) >= 0.001
            annotation('textbox',[mean([start_arrow add_arrow])-0.025 axes_pos(2)/3+axes_pos(4) 0.1 0.05],'String','**','FontSize',20,'FontWeight','demi','EdgeColor','none');
        elseif p(i) < 0.001
            annotation('textbox',[mean([start_arrow add_arrow])-0.03 axes_pos(2)/3+axes_pos(4) 0.1 0.05],'String','***','FontSize',20,'FontWeight','demi','EdgeColor','none');
        else
            annotation('textbox',[mean([start_arrow add_arrow])-0.03 axes_pos(2)/3+axes_pos(4) 0.2 0.05],...
                'String','ns','FontSize',14,'FontWeight','demi','EdgeColor','none');
%             annotation('textbox',[mean([start_arrow add_arrow])-0.07 axes_pos(2)/3+axes_pos(4) 0.2 0.05],'String',strcat('p =',sprintf(' %0.3f',p(i))),'FontSize',8,'FontWeight','demi','EdgeColor','none');
        end
        annotation('doublearrow',[start_arrow add_arrow],[axes_pos(2)/4+axes_pos(4) axes_pos(2)/4+axes_pos(4)]);

        start_arrow = start_arrow + (axes_pos(3)/(size(val,1)));
        add_arrow = add_arrow + (axes_pos(3)/(size(val,1)));
    end 
else
    row = 1;
    for i = 1:size(val_inp,2):size(val,1)
        if strcmp(usetest,'default')
            for j = i:i+size(val_inp,2)-2
                if (sum(~isnan(val(j,:)))>4 & sum(~isnan(val(j+1,:)))>4)
                    if ((lillietest(val(j,:))==0 & lillietest(val(j+1,:))==0 | ...
                                sum(~isnan(val(j,:)))>30 & sum(~isnan(val(j+1,:)))>30))
                        [~,p(row)] = ttest2(val(j,:),val(j+1,:),0.05,usetail);
                    else
                        p(row) = ranksum(val(j,:),val(j+1,:)); 
                    end
                else
                    if (sum(~isnan(val(j,:)))>3 & sum(~isnan(val(j+1,:)))>3)
                        p(row) = ranksum(val(j,:),val(j+1,:));
                    else
                        p(row) = NaN;
                    end
                end
                row = row + 1;
            end
        elseif strcmp(usetest,'param')
            for j = i:i+size(val_inp,2)-2
                [~,p(row)] = ttest2(val(j,:),val(j+1,:),0.05,usetail);
                row = row + 1;
            end
        elseif strcmp(usetest,'nonparam')
            for j = i:i+size(val_inp,2)-2
                if (sum(~isnan(val(j,:)))>3 & sum(~isnan(val(j+1,:)))>3)
                    p(row) = ranksum(val(j,:),val(j+1,:));
                else
                    p(row) = NaN;
                end
                row = row + 1;
            end
        end
    end
    
    for count = 1:length(p)
        if p(count) <= 0.05 & p(count) >= 0.01
            annotation('textbox',[mean([start_arrow add_arrow])-0.02 axes_pos(2)/3+axes_pos(4) 0.1 0.05],...
                'String','*','FontSize',20,'FontWeight','demi','EdgeColor','none');
%             annotation('textbox',[mean([start_arrow add_arrow])-0.07 axes_pos(2)/3+axes_pos(4) 0.2 0.05],...
%                 'String',strcat('p =',sprintf(' %0.3f',p(count))),'FontSize',10,'FontWeight','demi','EdgeColor','none');
        elseif p(count) < 0.01 & p(count) >= 0.001
            annotation('textbox',[mean([start_arrow add_arrow])-0.025 axes_pos(2)/3+axes_pos(4) 0.1 0.05],...
                'String','**','FontSize',20,'FontWeight','demi','EdgeColor','none');
%             annotation('textbox',[mean([start_arrow add_arrow])-0.07 axes_pos(2)/3+axes_pos(4) 0.2 0.05],...
%                 'String',strcat('p =',sprintf(' %0.3f',p(count))),'FontSize',10,'FontWeight','demi','EdgeColor','none');
        elseif p(count) < 0.001
            annotation('textbox',[mean([start_arrow add_arrow])-0.03 axes_pos(2)/3+axes_pos(4) 0.1 0.05],...
                'String','***','FontSize',20,'FontWeight','demi','EdgeColor','none');
        else
            annotation('textbox',[mean([start_arrow add_arrow])-0.025 axes_pos(2)/3+axes_pos(4) 0.2 0.05],...
                'String','ns','FontSize',14,'FontWeight','demi','EdgeColor','none'); % strcat('p=',sprintf('%0.3f',p(count)))
        end
        annotation('doublearrow',[start_arrow add_arrow],[axes_pos(2)/4+axes_pos(4) axes_pos(2)/4+axes_pos(4)]);
        
        if size(val_inp,2)>2
            if mod(count,size(val_inp,2)-1)==0
                start_arrow = start_arrow + 1.6/size(val,1);
                add_arrow = add_arrow + 1.6/size(val,1);
            else
                start_arrow = start_arrow + 0.8/size(val,1);
                add_arrow = add_arrow + 0.8/size(val,1);
            end
        else
            start_arrow = start_arrow + (axes_pos(3)/(size(val,1))) + 0.8/size(val,1);
            add_arrow = add_arrow + (axes_pos(3)/(size(val,1))) + 0.8/size(val,1);
        end
    end
    
end

% --- Resize y axis ---
ylimits = ylim;
ylim([nanmin(nanmin(val))-(ylimits(2)-ylimits(1))/10 nanmax(nanmax(val))+(ylimits(2)-ylimits(1))/10])

% --- Resize x axis ---
xlimits = xlim;
if plot_matrix == 0
    xrange = max(1-xlimits(1),xlimits(2)-size(val_inp,2));
else
    xrange = max(1-xlimits(1),xlimits(2)-size(val_inp,1));
end
if xrange<0.5
    xrange = 0.5;
end
if plot_matrix == 0
    xlim([1-xrange size(val_inp,2)+xrange])
else
    xlim([1-xrange size(val_inp,1)+xrange])
end

% --- Assign optional output ---
if nargout == 1
    varargout{1} = h;
end

end