function obj=stat_summary(obj,varargin)
% stat_summary Display summarized data for each value of X
%
% Example syntax (default arguments): gramm_object.stat_summary('type','ci','geom','lines','setylim',false)
% For each unique value of x, this can display various estimates of
% the location and variability of the corresponding y distribution.
% The optional 'name',value pairs can ne the following:
% - 'type':
%       - 'ci' : display the mean and the 95% confidence
%       interval of the mean (based on the assumption of a normal
%       distribution
%       - 'bootci' : display the mean and the 95% confidence
%       interval of the mean computed by bootstrap
%       - 'sem' : display the mean and the standard error of
%       the mean
%       - 'std' : display the mean and the standard deviation
%       and the mean
%       - 'quartile': display the 25% 50% (median) and 75%
%       percentiles
%       - '95percentile': display the median and 2.5 and 97.5
%       percentiles
%        - 'fitnormalci'
%       - 'fitpoissonci'
%        - 'fit95percentile'
% - 'geom':
%       - 'line': displays a line that connects the central locations
%       (mean,median)
%       - 'lines': displays a line that connects the central locations
%       and lighter lines for the variabilities
%       - 'area': displays a line that connects the central locations
%       and a transparent area for the variabilities. WARNING:
%       this changes the renderer to opengl and disables proper
%       vector output on older matlab versions
%       - 'solid_area': displays a line that connects the locations
%       and a solid area for the variabilities. Use this for
%       export to vector output.
%       - 'errorbar': displays error bars for variabilities.
%       - 'black_errorbar': displays black error bars for variabilities.
%       - 'bar': displays the locations as bars
% - 'setylim': set to true if you want the y axis limits to be
% set by the summarized data instead of the underlying data
% points.
% - 'interp': Use to interpolate the output, takes the same parameters as
% interp1 in order to specify the interpolation type. When the
% polar mode is specified as closed, or when 'interp','polar' is used
% the interpolation uses interpft, which supposes regular sampling around the circle.
% - 'interp_in': Use to (linearly) interpolate the input. This is intended
% for input given as cells when the x value is different for
% each cell and not aligned. The argument corresponds to the
% number of points used to generate the interpolation. Ideally
% for this the number of number of points should be higher than
% the x resolution of the data, otherwise some data will be
% unused.
% - 'bin_in': Use to bin the input. This is intended for input
% given as a 1-D array, creates bins over x and computes the summary
% over the binned data. Argument corresponds to the number of
% bins
% - 'dodge': use to dodge the plotted elements depending on
% color (recommended for 'bar', 'errorbar', 'black_errorbar').
% A value of 0 deactivates dodging. Other values set the space
% between the dodged elements as ratio of the x intervals.
% - 'width': use to set the width of bars and error bars (error
% bars are 1/4 th the width of bars).
%
% A setting of 'dodge',1,'width',1 will create bars that are
% fully dodged (ie don't overlap, but are not separated) and where
% all bars occupy the full interval between the x values, e.g.:
% * __    *
% *|  |__ *
% *|  |  |*
% *|__|__|*
%
% A setting of 'dodge',0.8,'width',0.8 will have fully dodged bars that are
% not separated, but only occupy 50% of the space between x
% values, e.g.:
% *  _    *
% * | |_  *
% * | | | *
% * |_|_| *
%
%A setting of 'dodge',0.8,'width',0.6 will add some
% spacing between the bars, e.g.:
% *  _      *
% * | |  _  *
% * | | | | *
% * |_| |_| *
%
% A setting of 'dodge',0.8,'width',1 will create dodged but overlapping
% bars, e.g.:
% *  ___    *
% * |  _|_  *
% * | |   | *
% * |_|___| *



p=inputParser;
my_addParameter(p,'type','ci'); %'95percentile'
my_addParameter(p,'geom','area');
my_addParameter(p,'dodge',[]);
my_addParameter(p,'width',[]);
my_addParameter(p,'setylim',false);
my_addParameter(p,'interp','none');
my_addParameter(p,'interp_in',-1);
my_addParameter(p,'bin_in',-1);
parse(p,varargin{:});

obj.geom=vertcat(obj.geom,{@(dd)my_summary(obj,dd,p.Results)});
obj.results.stat_summary={};
end




function hndl=my_summary(obj,draw_data,params)

%Advanced defaults
if isempty(params.dodge)
    if sum(strcmp(params.geom,'bar'))>0 && draw_data.n_colors>1 %If we have a bar as geom, we dodge
        params.dodge=0.6;
    else
        params.dodge=0;
    end
end

if isempty(params.width) %If no width given
    if params.dodge>0 %Equal to dodge if dodge given
        params.width=params.dodge*0.8;
    else
        params.width=0.5;
    end
end



if iscell(draw_data.x) || iscell(draw_data.y) %If input was provided as cell/matrix
    
    if params.interp_in>0
        %If requested we interpolate the input
        uni_x=linspace(obj.var_lim.minx,obj.var_lim.maxx,params.interp_in);
        [x,y]=cellfun(@(x,y)deal(uni_x,interp1(x,y,uni_x,'linear')),draw_data.x,draw_data.y,'UniformOutput',false);
        y=padded_cell2mat(y);
    else
        %If not we just make a padded matrix for fast
        %computations (we'll assume that X are roughly at the
        %same location for the same indices)
        x=padded_cell2mat(draw_data.x);
        y=padded_cell2mat(draw_data.y);
        uni_x=nanmean(x);
    end
    
    if params.bin_in>0
        warning('bin_in in stat_summary() not supported for Matrix/Cell X/Y inputs');
    end
    
    if strfind(params.type,'fit')
        %If we have a params.type using distributions fits we
        %can't vectorize the call to computeci so we do it in a for
        %loop
        ymean=zeros(length(uni_x),1);
        yci=zeros(length(uni_x),2);
        for ind_x=1:length(uni_x)
            [ymean(ind_x),yci(ind_x,:)]=computeci(y(:,ind_x),params.type);
        end
    else
        [ymean,yci]=computeci(y,params.type);
    end
    
else %If input was provided as 1D array
    
    x=comb(draw_data.x);
    y=comb(draw_data.y);
    
    if params.bin_in>0
        %If X binning was requested we do it
        binranges=linspace(obj.var_lim.minx,obj.var_lim.maxx,params.bin_in+1);
        bincenters=(binranges(1:(end-1))+binranges(2:end))/2;
        [~,binind]=my_histcounts(x,binranges,'count');
        uni_x=bincenters;
        sel=binind~=0; %histcounts can return zero as bin index if NaN data we remove them here
        x=bincenters(binind(sel));
        y=y(sel);
    else
        if sum(strcmp(params.geom,'area'))>0 || sum(strcmp(params.geom,'line'))>0 ||...
                sum(strcmp(params.geom,'lines'))>0  || sum(strcmp(params.geom,'solid_area'))>0
            %To avoid interruptions in line and area plots we
            %compute uniques over current data only
            uni_x=unique(x); %Sorted is the default
        else
            %compute unique Xs at the facet level to avoid
            %weird bar sizing issues when dodging and when
            %colors are missing
            facet_x=comb(draw_data.facet_x);
            uni_x=unique(facet_x);
        end
        %Here we need to implement a loose 'unique' because of
        %potential numerical errors
        uni_x(diff(uni_x)<1e-10)=[];
    end
    
    if params.interp_in>0
        warning('inter_in in stat_summary() not supported for non Matrix/Cell X/Y inputs');
    end
    
    ymean=nan(length(uni_x),1);
    yci=nan(length(uni_x),2);
    
    %Loop over unique X values
    for ind_x=1:length(uni_x)
        %And here we have a loose selection also because of
        %potential numerical errors
        ysel=y(abs(x-uni_x(ind_x))<1e-10);
        
        if ~isempty(ysel)
            [ymean(ind_x),yci(ind_x,:)]=computeci(ysel,params.type);
        end
    end
end


%Do we set the y limits according to the smoothed curves or to
%the original data ?
if params.setylim
    if sum(sum(isnan(yci)))~=numel(yci) %We only do this if yci is not weird
        if obj.firstrun(obj.current_row,obj.current_column) %Initialize for the first run in the subplot
            obj.plot_lim.maxy(obj.current_row,obj.current_column)=max(max(yci));
            obj.plot_lim.miny(obj.current_row,obj.current_column)=min(min(yci));
        else %Update for subsequent runs in the subplot
            if max(max(yci))>obj.plot_lim.maxy(obj.current_row,obj.current_column)
                obj.plot_lim.maxy(obj.current_row,obj.current_column)=max(max(yci));
            end
            if min(min(yci))<obj.plot_lim.miny(obj.current_row,obj.current_column)
                obj.plot_lim.miny(obj.current_row,obj.current_column)=min(min(yci));
            end
        end
    end
end

%When we do bar plots we want to have zero in the y axis anyway
if sum(strcmp(params.geom,'bar'))>0
    if obj.plot_lim.miny(obj.current_row,obj.current_column)>0 %Values above zero -> change miny
        obj.plot_lim.miny(obj.current_row,obj.current_column)=0;
    end
    if obj.plot_lim.maxy(obj.current_row,obj.current_column)<0 %Values below zero -> change maxy
        obj.plot_lim.maxy(obj.current_row,obj.current_column)=0;
    end
end

%Do we interpolate the summary results for display ?
if ~strcmp(params.interp,'none')
    if size(yci,1)>2
        yci=yci';
    end
    
    if obj.polar.is_polar && obj.polar.is_polar_closed && ~strcmp(params.interp,'polar')
        disp([params.interp ' interpolation overriden, ''polar'' used']);
        params.interp='polar'; %%If the plot is polar we override the interpolation type do an optimal fft interpolation
    end
    
    if strcmp(params.interp,'polar')
        uni_x=0:pi/50:2*pi-pi/50;
        ymean=interpft(ymean,100);
        tmp_yci1=interpft(yci(1,:),100);
        tmp_yci2=interpft(yci(2,:),100);
        yci=[tmp_yci1 ; tmp_yci2];
    else
        %For non polar plots we do a regular interpolation
        new_x=linspace(min(uni_x),max(uni_x),100);
        ymean=interp1(uni_x,ymean,new_x,params.interp);
        tmp_yci1=interp1(uni_x,yci(1,:),new_x,params.interp);
        tmp_yci2=interp1(uni_x,yci(2,:),new_x,params.interp);
        yci=[tmp_yci1 ; tmp_yci2];
        uni_x=new_x;
    end
end

%Store results
obj.results.stat_summary{obj.result_ind,1}.x=uni_x;
obj.results.stat_summary{obj.result_ind,1}.y=ymean;
obj.results.stat_summary{obj.result_ind,1}.yci=yci;

%Do the actual plotting
hndl=plotci(obj,uni_x,ymean,yci,draw_data,params.geom,params.dodge,params.width);

%Copy handles
hnames=fieldnames(hndl);
for k=1:length(hnames)
    obj.results.stat_summary{obj.result_ind,1}.(hnames{k})=hndl.(hnames{k});
end
end



function [ymean,yci]=computeci(y,type)

ymean=nanmean(y);
try
    switch type
        case 'bootci'
            yci=bootci(200,@(y)nanmean(y),y);
        case 'ci'
            ci=1.96*nanstd(y)./sqrt(sum(~isnan(y)));
            yci=bsxfun(@plus,ymean,[-ci;ci]);
        case 'std'
            ci=nanstd(y);
            yci=bsxfun(@plus,ymean,[-ci;ci]);
        case 'sem'
            ci=nanstd(y)./sqrt(sum(~isnan(y)));
            yci=bsxfun(@plus,ymean,[-ci;ci]);
        case 'quartile'
            ymean=nanmedian(y);
            yci=prctile(y,[25 75]);
        case '95percentile'
            ymean=nanmedian(y);
            yci=prctile(y,[2.5 97.5]);
        case 'fitnormalci'
            pd=fitdist(y,'Normal');
            ymean=pd.mean();
            ci=pd.paramci();
            yci=ci(:,1)';
        case 'fitpoissonci'
            pd=fitdist(y,'Poisson');
            ymean=pd.mean();
            ci=pd.paramci();
            yci=ci(:,1)';
        case 'fitbinomialci'
            pd=fitdist(y,'Binomial');
            ymean=pd.mean;
            ci=pd.paramci;
            yci=ci(:,2)';
        case 'fit95percentile'
            pd=fitdist(y,'Normal');
            ymean=pd.icdf(0.5);
            yci=pd.icdf([0.025 0.975]);
        otherwise
            warning(['Unknown CI type ' type]);
    end
catch
    disp('Not enough samples for CI computation...skipping')
    yci=repmat([NaN NaN],size(y,2));
end
end